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ABSTRACT 


This  thesis  applies  time-frequency  transformations  to  radar  signals.  Specifically,  it 
considers  the  feasibility  of  applying  time-frequency  transformations  to  extract  the  intra- 
pulse  modulation  parameters  of  radar  signals.  In  this  work,  we  consider  radar  signals  with 
analog  pulse  compression;  specifically  linear  or  hyperbolic  intra-pulse  modulation. 
Several  time-frequency  transformations  are  investigated  to  identify  which  one  gives  the 
most  accurate  image  representation  for  signals  in  noisy  environments.  Next,  image 
processing  techniques  are  applied  in  conjunction  with  an  adaptive  curve  fitting  method, 
for  the  hyperbolic  modulation  scheme,  to  extract  the  parameters  of  the  frequency 
equation.  Results  show  that  for  the  linear  chip  case  the  frequency  equation  can  be 
estimated  with  small  error  down  to  SNR  equal  to  — lOdB.  The  proposed  method  for  the 
hyperbolic  chirp  modulation  is  less  immune  to  noise  degradation  and  it  can  be  used  down 
to  SNR  level  equal  to  2dB. 
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I.  INTRODUCTION 


The  analysis  of  a  signal  in  both  the  time  and  frequency  domain  gives  a  more 
complete  description  than  that  provided  by  individual  analysis  in  either  domain.  This 
advantage  was  recognized  a  long  time  ago  but  the  complexity  of  the  joint  time-frequency 
algorithms  made  their  implementations  an  extremely  slow  operation  due  to  computing 
power  limitations.  However,  as  computers  become  more  powerful,  joint  time-frequency 
analysis  may  be  applied  in  a  larger  number  of  real-world  problems. 

A.  OBJECTIVE 

This  thesis  investigates  the  application  of  time-frequency  transformations  to  the 
extraction  of  intra-pulse  modulation  parameters  from  radar  signals.  This  information  can 
be  a  useful  tool  to  identify  the  specific  radar  type  and  help  with  jamming  techniques,  if  so 
desired.  Intra-pulse  modulation  parameters  are  usually  extracted  using  hardware  schemes 
which  are  well  suited  to  extract  the  instantaneous  frequency  at  high  SNR  levels. 

Our  study  considers  the  problem  from  a  different  angle  and  investigates  the 
application  of  time-frequency  and  a  basic  image  processing  technique  to  extract  the 
information.  Joint  time-frequency  transformations  can  be  implemented  inexpensively  and 
give  accurate  results  in  medium  to  high  SNR  levels.  We  consider  radar  signals  with 
analog  pulse  compression;  specifically  linear  or  hyperbolic  intra-pulse  modulation. 
Several  time-frequency  transformations  are  investigated  to  identify  that  which  gives  the 
most  accurate  image  representation  for  signals  in  noisy  environments.  Next,  image 
profcessing  techniques  are  applied  used  in  conjunction  with  an  adaptive  curve  fitting 
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method,  to  extract  the  parameters  of  the  frequency  equation.  The  type  of  modulation,  as 
well  as,  the  start  and  stop  time  of  the  pulse  under  investigation  was  assumed  to  be  known. 

B.  THESIS  ORGANIZATION 

Radar  signals  with  large  time  bandwidth  products  are  introduced  in  chapter  H. 
Chapter  IH  presents  an  overview  of  the  time  frequency  methods  considered.  Chapter  IV 
discusses  the  basic  idea  behind  the  Radon  transform  and  its  application  in  detecting  line 
parameters  from  a  noisy  image.  Chapters  V  and  VI  present  the  complete  methods  and  the 
simulaitions  conducted  to  test  the  schemes  derived  to  extract  the  intra-pulse  modulation 
parameters  for  linear  and  hyperbolic  modulation.  Finally,  conclusions  and 
recommendations  for  further  research  are  presented  in  chapter  VII. 
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11.  RADAR  SIGNALS  WITH  LARGE  TIME  BANWIDTH  PRODUCTS 


Ideally  we  would  like  a  radar  to  produce  good  range  resolution  as  well  as  large 
detection  range.  Unfortunately,  these  requirements  are  contradictory  since  the  first 
requires  a  small  pulse  period,  while  the  second  requires  a  large  pulse  duration.  The 
problem  can  be  solved  with  a  technique  that  uses  intra-pulse  modulation  or,  as  it  is  better 
known,  pulse  compression.  This  chapter  present  the  basic  ideas  behind  pulse 
compression. 


A.  RESOLUTION 

A  very  important  aspect  of  the  radar  is  its  resolution  properties.  By  this  we  mean 
the  ability  to  detect  multiple  targets  that  are  close  together.  Generally,  there  are  four  types 
of  resolution: 

-  Range  resolution 

-  Horizontal  (azimuth)  cross-range  resolution 

-  Vertical  (elevation)  cross-range  resolution 

-  Doppler  frequency  resolution. 


Figure  1:  Range  and  Cross-Range  Resolution  From  [1] 
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1.  Range  Resolution 

Range  resolution  is  the  ability  to  distinguish  between  targets,  which  are  at  the 
same  angular  position  but  are  at  different  distances  to  the  radar.  This  resolution  is  directly 
proportional  to  the  radar  signal  bandwidth:  the  larger  the  bandwidth  is,  the  better  the 
range  becomes.  In  order  to  resolve  two  targets,  the  basic  criterion  is  that  they  must  be 
separated  by  at  least  the  range  equivalent  of  the  width  of  the  processed  echo  pulse  as 
illustrated  in  Figure  2. 


Processed  Unresolved  Near  resolved  Resolved 

pulse  width 


a.  Wide  processed  pulse  width 


Processed  Unresolved  Near  resolved  Resolved 

pulse  width 


Figure  2:  Range  Resolution  and  Processed  Pulse  Width.  From  Ref.  [1 J 

2.  Horizontal  and  Vertical  Cross-Range  Resolution. 

Cross-range  resolution  is  defined  as  the  ability  to  distinguish  between  targets  that 
have  the  same  distance  from  the  radar  but  are  at  a  different  azimuth  (horizontal)  and/or 
different  elevation  (vertical).  This  resolution  is  inversely  proportional  to  the  antenna 
beamwidth.  That  means  the  smaller  the  beamwidth  is,  the  better  the  cross-range 
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resolution  is.  The  basic  constraint  is  that  the  targets  must  be  separated  by  at  least  the 


antenna  beamwidth.  Figure  3  shows  the  effect  of  the  antenna’s  beamwidth  to  the 
horizontal  cross-range  resolution. 


Antenna 


scan  _ 

A  k  _ _ ^ 

F*^  A  k  ^ 

^•4  ^ 

Targets  separoted  by 
less  than  beamwidth 


Targets  separated  by 
antenna  beamwidth 


Targets  separated  by 
greater  than  beamwidth 


Signal  echoes 
Unresolved 


Signal  echoes  - 
Near  resolved 


Signal  echoes 
Resolved 


Figure  3:  Antenna  Beamwidth  and  Horizontal  Cross-Range  Resolution.  From 
Ref.  [1] 


3.  Doppler  Resolution 

Radar  doppler  resolution  is  defined  as  the  ability  to  distinguish  between  targets 
that  have  the  same  range,  azimuth  and  elevation  but  have  different  radial  velocities.  The 
basic  criterion  is  that  the  Doppler  frequencies  of  the  targets  must  be  separated  by  at  least 
one  cycle  over  the  time  of  observation. 


B.  PULSE  COMPRESSION 

A  radar  should  have  a  large  a  detection  range.  If  we  assume  that  we  use  the  best 
techniques  for  antenna  and  signal  processing  gain,  the  only  way  to  increase  the  range  is 
by  increasing  the  energy  of  the  transmitted  signal.  There  are  two  ways  to  increase  the 
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energy;  by  increasing  its  amplitude  (power)  of  the  pulse,  or  by  increasing  the  pulse 
duration.  Since  all  transmitters  have  peak  power  limitations,  the  most  common  choice  is 
to  increase  the  pulse  duration,  which  is  in  contradiction  to  the  need  for  narrow  pulses  to 
ensure  good  range  resolution.  The  technique  that  uses  long  pulses  with  high  bandwidth  is 
called  pulse  compression. 

1.  Analog  Pulse  Compression 

This  technique  uses  modulated  radio  frequency  signals.  Let’s  assume  that  we  have 
a  finite  pulse  duration  radar  signal  with  a  constant  radio  frequency  defined  as: 

c(r)  =  r(0*5(0,  (H-l) 

where  r(t)  is  a  pulse  and  s(t)  is  a  sinusoid  function.  Let  C(f),  R(f)  and  S(f)  be  the 
Fourier  transform  of  the  time  domain  signals,  then: 

C(f)  =  R(f)*S(f),  01-2) 

where  the  operator  (*)  denotes  the  convolution  operation.  Thus,  the  spectrum  of 
the  total  signal  has  the  form  of  a  sine  function  centered  at  the  frequency  of  the  sinusoid 
with  a  null  to  null  bandwidth  equ^  to  the  double  of  the  reciprocal  of  the  width  of  the 
pulse  in  the  time  domain,  as  illustrated  in  Figure  4. 
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Figure  4:  Bandwidth  of  a  CW  pulse. 


The  Fourier  transform  S(f)  is  not  a  delta  function  when  the  radio  signal  5(0  is  a 
modulated  sinusoid.  This  in  turns  affects  the  bandwidth  of  the  transmitted  signal  and 
gives  a  null  to  null  bandwidth  greater  than  2/T,  as  shown  in  Figure  5. 


Figure  5:  Bandwidth  of  a  linear  modulated  pulse. 
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Next  we  will  examine  the  two  most  widely  used  types  of  pulse  intra-modulation; 
linear  and  hyperbolic  modulations. 

a)  Linear  Modulation 

In  this  case  the  frequency  of  the  radio  signal  inside  the  pulse  changes 
linearly  with  a  slope  k  from  an  initial  frequency  fo  which  leads  to; 

f(t)  =  fo+k-t.  (n-3) 

Where /o  is  the  initial  frequency  at  t=0.Thus,  the  transmitted  signal  c(t)  is 

given  as: 

c(t)  =  rect  C-^)  •  exp(  y 2;r J  fdt )  =rect  C-^)  •  >  .  (n-4) 

b)  Hyperbolic  Modulation 

In  this  case,  the  signal  inside  the  pulse  is  modulated  following  a  hyperbolic 
equation.  The  resulting  frequency  equation  is  given  as: 

f{t)=j  +  t>.  (n-5) 

Note  that  the  above  time  varying  frequency  f(t)  is  not  finite  for  t=0.  Thus, 
the  above  definition  is  modified  by  introducing  a  time  shift  to  to  ensure  a  finite  frequency 
value  for  t=0,  which  leads  to: 

/(0  =  -^+^>.  (n-6) 

l+tQ 
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Thus,  the  instantaneous  frequency  yftj  can  be  viewed  as  being  a  part  of  the 


hyperbola  defined  from  the  equation  (11-5)  with  an  initial  value  of  = 


Time 

Figure  6:  Hyperbolic  frequency  modulation. 


As  a  result,  the  equation  of  the  transmitted  signal  is  given  by: 
c(r)  =  rect (j)  •  exp(  jlTrj  fdt)  =rect  (^)  • 


2.  Digital  Pulse  Compression 

Digital  pulse  compression  involves  phase-coded  waveforms.  These  are  usually  bi¬ 
phase  modulated  sinusoids  with  the  two  possible  phases  to  be  at  0°  and  180°.  The  overall 
digital  waveform  consists  of  an  array  of  N  subpulses,  each  one  with  an  initial  phase  of  0° 
or  180°  and  a  time  length  Ts.  The  width  of  the  total  waveform  is  Te-  Figure  7  shows  a 
biphase  coded  signal  using  a  13-bit  biphase  sequence.  In  this  figure  a  “+”  corresponds  to 
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a  subpulse  with  an  initial  phase  0°  and  a  corresponds  to  a  subpulse  with  an  initial 
phase  180°. 


a.  Waveform 


± 

+  +  + 

S 

□ 

B 

l+l+i 

-  +  - 

FI 

- 

•^s 

Tc. 

- - - 

Figure  7:  Phase  coded  waveform.  From  Ref.fl] 

The  codes  that  are  used  in  the  phase  coded  waveforms  should  have  a  delta-like 
autocorrelation  function,  to  allow  pulse  compression  to  take  place  without  windowing  the 
signal.  The  two  most  widely  used  codes  are  the  Barker  codes  and  the  Pseudorandom 
codes  [1]. 
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III.  TIME-FREQUENCY  REPRESENTATIONS 

Generally,  a  signal  can  be  represented  using  different  types  of  decompositions. 
The  two  most  widely  used  representations  are  the  time  domain  and  the  frequency  domain 
representation.  The  first  shows  how  the  amplitude  of  the  signal  changes  with  respect  to 
time  while  the  second  shows  how  often  these  changes  occur.  These  representations  are 
uniquely  related  to  each  other  via  the  Fourier  transform. 

A.  FOURIER  ANALYSIS 

1.  Fourier  Series 

The  basic  idea  behind  the  Fourier  decomposition  is  to  express  a  periodic  signal  as 
a  weighted  sum  of  sinusoids.  Let’s  assume  that  we  have  a  signal  x(t)  with  period  T. 
Provided  that  it  satisfies  the  Dirichlet’s  conditions,  then  it  can  represented  as  follows: 

(m-i) 

*=-~ 

where 

,  =— .•  (in-2) 

^  T 

The  frequency  1/T  is  called  the  fundamental  frequency.  The  equation  (III-l)  gives 
the  general  case  where  the  type  of  basis  functions  are  complex  exponentials. 

2.  Fourier  Transform 

By  analogy  with  equation  (111-2),  the  Fourier  transform  X(f)  of  a  continuous  signal 
x(t)  is  given  by: 
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(ni-3) 


^(/)=  \x{t)-e-^^’^dt. 

—  oo 

The  plot  of  the  magnitude  squared  of  the  Fourier  transforai  is  proportional  to  the 
spectrum  of  the  signal  x(t).  The  original  signal  can  be  recovered  from  its  Fourier 
transform  with  the  inverse  Fourier  transform: 

oo 

x{t)=  .  (in-4) 

—  oo 

The  fast  discrete  version  of  Fourier  transform  revolutionized  the  signal  processing 
area.  However,  there  exist  many  signals  which  cannot  be  described  accurately  with  this 
transformation.  Although  the  spectrum  indicates  the  frequencies  and  intensities  of  the 
signal,  it  provides  no  information  regarding  the  variations  of  the  signal  characteristics  as  a 
function  of  time.  Thus,  the  Fourier  transform  does  not  analyze  non  stationary  signals  very 
accurately.  For  example,  consider  two  signals  Xi(t)  and  X2(t),  observed  during  the  same 
time  period  T,  as  shown  in  Figure  8.  The  first  one  is  the  sum  of  two  pure  tones  while  the 
second  one  consists  of  the  first  tone  for  the  first  half  time  and  then  switches  to  the  second 
tone  for  the  rest  of  its  duration.  Their  Fourier  transform  magnitudes  look  almost  identical, 
however  they  are  not  exactly  the  same  since  the  Fourier  transform  is  a  reversible 
operation.  Figure  8  shows  that  it  is  not  possible  to  extract  specific  information  regarding 
the  frequency  variations  as  a  function  of  time.  Unfortunately,  most  signals  in  real  life  are 
not  stationary,  such  as  for  example  speech,  music,  vibration  signals  and  many  more. 
Analyzing  non-stationary  signals  requires  a  tool  which  allows  to  represent  the  variations 
of  the  frequency  content  as  a  function  of  time;  i.e.,  a  joint  time-frequency  representation. 
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Several  techniques  have  been  developed  for  this  type  of  representation.  In  this  thesis,  we 
will  deal  with  two  major  classes:  atomic  decompositions  and  energy  distributions. 


x1(t) 


Figure  8:  Fourier  transform  of  two  signals. 

B.  ATOMIC  DECOMPOSITIONS 

The  Fourier  transform  can  be  considered  as  the  projection  of  a  signal  into  an 
infinite  set  of  sinusoids.  The  problem  with  sinusoids  is  that  they  are  not  localized  in  time, 
although  they  are  perfectly  localized  in  frequency.  The  basic  idea  behind  the  atomic 
decomposition  is  to  decompose  a  signal  into  a  set  of  atoms  (i.e.  functions)  that  are  well 
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localized  both  in  time  and  in  frequency.  These  types  of  decompositions  are  called  atomic 
decompositions  or  linear  time-frequency  representations. 

1.  Short-Time  Fourier  Transform 

This  method  is  an  extension  of  the  Fourier  transform,  where  localization  in  time  is 
introduced  by  windowing  the  sinusoidal  functions  used  in  the  decomposition.  Thus,  the 
Sort-Time  Fourier  Transform  (STFT)  is  given  by: 

OO 

(in-5) 

— OO 

where  g{t)  is  the  window  function.  Note  that  the  STFT  is  an  invertible  operation, 
provided  that  the  window  is  of  finite  energy,  i.e., 

.|  OO  OO 

xit)  =  —  •  J  J S{T,  f)  ■  git  -  T)  ■  •  df  ,  (m-6) 

-00-00 

OO 

where  Ef^  =  The  STFT  operation  can  be  considered  as  a  set  of  successive 

— OO 

Fourier  transforms  applied  to  windowed  segments  of  the  signal,  as  illustrated  in  Figure  9. 
Thus,  g(t)  can  be  viewed  as  a  sliding  window  which  allows  for  segmentation  of  the 
original  signal. 
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Figure  9:  Short-Time  Fourier  Transform.  From  Ref  [2] 

The  STFT  of  x(t)  can  also  be  seen  as  the  expansion  of  the  signal  into  a  set  of 
atomic  functions  of  the  form: 

=  (in-7) 

The  above  function  can  be  recognized  as  a  windowed  complex  sinusoid,  where 
the  window  type  and  duration  affects  the  resolution  of  the  STFT  in  time  and  in  frequency. 
The  plot  of  the  squared  magnitude  of  the  STFT  is  called  the  spectrogram. 

2.  Gabor  Expansion 

In  1946,  Dennis  Gabor  introduced  the  decomposition  of  a  signal  into  a  weighted 
sum  of  functions  that  are  both  localized  in  time  and  frequency.  The  Gabor  expansion  is 
defined  as: 

oo  oo 

(in-8) 
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where: 


it)  =  hit -mT)-  exp(  jnQ.t) .  (111-9) 

The  coefficients  Cni,n  CElled  the  Gabor  coefficients,  while  T  and  ^2  are  the  time 
and  frequency  sampling  steps.  Initially,  Gabor  selected  h(t)  to  be  Gaussian,  since  it  is 
optimally  concentrated  in  the  joint  time-frequency  domain.  However,  h(t)  can  be  any 
function.  Comparing  equations  (III-8)  and  (III-6)  shows  that  the  Gabor  expansion  is 
similar  to  the  discrete  version  of  the  inverse  STFT.  In  fact,  the  Gabor  expansion  is  the 
generalization  of  the  discrete  inverse  STFT  which  is  used  as  a  fast  algorithm  for  its 
computation  [2,9]. 

3.  Wavelet  Analysis 

a)  Continuous  wavelet  transform 

We  saw  that  the  STFT  of  x(t)  is  the  decomposition  of  the  signal  into  a  set 
of  atoms  of  the  form  (III-7)  which  are  two-parameter  functions  of  a  complex  sinusoid. 
The  continuous  wavelet  transform  is  a  more  general  case,  where  the  basis  function  is  not 
a  complex  sinusoid  but  a  function  with  specific  properties  of  the  form: 

=  (m-lO) 

The  parameters  a  and  r  are  the  scale  and  translation  factor,  respectively. 
The  function  'P(t)  is  called  the  mother  wavelet.  This  basis  function  must  meet  two 
important  criteria;  1)  it  must  have  finite  duration,  and  2)  it  must  have  a  zero  average 
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value.  Some  of  the  most  widely  used  wavelets  are  the  Haar,  Daubechies,  Coiflet, 
Symmlet  which  are  shown  in  Figure  10. 


Haar 


Daubechie-4 


Figure  10:  Different  types  of  wavelets.  From  [12] 


The  continuous  wavelet  transform  (CWT)  of  x(t)  is  defined  as  the  integral 
of  the  signal  multiplied  by  scaled  and  shifted  versions  of  the  mother  wavelet  function  \|r(t) 
[6]: 


C(T,a)  =  ^-]xit)-W\^)dt , 


(m-ll) 


where  the  factor  normalizes  the  transformation.  Comparing  equation  (III-ll)  with 

equation  (III-5)  shows  that  we  can  relate  the  scale  (or  dilation)  factor  a  to  the  frequency/. 
The  difference  in  the  CWT  is  that  the  time  and  frequency  resolution  are  also  controlled 
by  the  factor  a  instead  of  the  window  function  only,  as  is  the  case  for  the  SIFT'. 
Specifically,  small  values  of  a  mean  that  the  wavelet  is  contracted  in  time,  thus  expanded 
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in  frequency.  Therefore,  the  time  resolution  is  good  while  the  frequency  resolution  is 
poor.  Similarly  the  wavelet  is  expanded  in  time  and  contracted  in  frequency  when  a  takes 
large  values.  In  this  case,  the  time  resolution  becomes  poor  as  the  frequency  resolution 
improves.  This  unique  characteristic  of  the  wavelet  analysis  is  known  as  the 
multiresolution  capability  in  the  time-frequency  plane.  Figure  1 1  shows  the  tiling  of  the 
time-frequency  plane  for  the  wavelet  analysis  and  for  the  STFT. 


0  Time  -j 


Figure  11:  STFT  vs.  Wavelet  time-frequency  resolution. 
Equation  (IH-I  1)  can  also  be  written  in  the  form: 


C{a,r)  =  x{t)*^  a{t,r) , 

where  *  denotes  convolution  and 


(ffl-12) 


Va  Cl 


(in-13) 
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It  can  be  shown  that  'Pa(t)  is  the  impulse  response  of  a  band-pass  filter 
[4].  Thus,  equation  (111-12)  shows  that  the  CWT  of  x(t)  can  be  computed  by  filtering  the 
signal  through  a  series  of  band-pass  proportional  filters.  The  plot  of  the  square  magnitude 
of  C(a,  t)  is  called  the  scalogram,  by  analogy  with  the  STFT  spectrogram.  The  CWT  is  a 
reversible  operation,  which  means  that  the  original  signal  x(t)  can  be  derived  from  its 
transform  C(a,  t)  by  [4]: 


x(r)  =  — f  {  C(T,a)-^-=-y/(^—^)drda, 


(in-14) 


where 


CD 


dco  , 


(in-15) 


and  W(cd)  is  the  Fourier  transform  of  TCr) .  Equation  (111-14)  implies  that  C(a,T)  needs  to 
be  known  for  the  whole  range  of  a  in  [0,<»] ,  in  order  to  recover  the  signal  in  time  domain 
from  the  CWT.  When  the  CWT  is  known  for  values  a<ao  only,  a  scaling  function  which 
contains  the  information  for  the  range  a  >  needs  to  be  introduced.  The  scaling  function 

<t)(t)  is  the  impulse  response  of  a  lowpass  filter.  If  we  define: 

<z),(0=^-«>(-),  (in-16) 

•4s  s 
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then,  equation  (111-14)  can  be  written  as: 


1  °  d  \ 

X(t)  =  —  lc(a,r)*%(0-^+ - L(ao,T)*^^  (0,  (0-17) 

^yf  0  ^  ’^0  ° 

where  the  term  L(a,r)  is  called  the  low  frequency  approximation  of  x(t)  at  scale  a  and  is 
given  by: 


(m-18) 


b)  Discrete  wavelet  transform 

The  discrete  wavelet  transform  (DWT)  can  be  derived  from  the  continuous 
wavelet  transform  expression  using  the  discrete  versions  of  a,T  and  t.  Specifically,  the 
DWT  of  a  discrete  signal  x(n)  is  defined  as: 

C(.a,b)  =  ^-^x(.n)'v\—).  (ni-19) 

In  practice  we  further  restrict  the  factors  a  and  b  to 

a  =  2f  b  =  k-2^  ,  (m-20) 


where  k  and  j  are  integers.  Using  equation  (111-20),  equation  (III-19)  can  be  written  as: 

N 


CjX  =f,-frX(,ny¥\2-‘n-k). 

n=i  y2'^ 


(m-21) 


Equation  (111-21)  shows  that  the  DWT  is  decimated  by  a  factor  of  two  at 
each  successive  scale  j.  Specifically,  the  maximal  values  of  j  and  k  are  log2(N),  and  N2'^ 
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for  a  given  value  of  k,  for  a  signal  length  equal  to  N.  The  signal  x(n)  can  be  recovered 
from  the  DWT  using  the  discrete  version  of  equation  (HI- 17)  which  yields  [4]: 


In  2  A 

4n]  =  — + - Lj^j,  (10-22) 

j=l 

where  the  operator  0  denotes  the  circular  convolution.  Equation  (111-22)  shows  that  the 
signal  can  be  decomposed  in  the  coefficients  Cj,k  and  Ljo,k  using  a  pair  of  low-pass  and 
high-pass  filters.  This  scheme  was  first  proposed  by  Mallat  and  uses  the  quadrature 
mirror  filters  (QMF)  theory  [4].  Specifically,  the  signal  x(n)  is  fed  into  a  pair  of  low-pass 
(LP)  and  high-pass  (HP)  filters  which  cover  the  frequency  range  from  0  to/^ /2,  as  shown 
in  Figure  12.  The  LP  filter  output  contains  the  approximation  of  the  signal  and 
corresponds  to  the  first  term  of  equation  (111-22),  while  the  output  of  the  HP  filter  gives 
the  details  of  the  signal  and  corresponds  to  the  second  term  of  equation  (111-22).  The  filter 
outputs  can  be  decimated  by  a  factor  of  two,  since  each  filter  covers  only  half  of  the 
original  frequency  range. 


x(n) 


LP 


HP 


\  /■ 
\  ✓ 
A 


t/4 


m 


Figure  12:  Quadrature  mirror  filters 
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This  scheme  continues  at  the  next  level  by  decomposing  the 
approximation  into  second  stage  details  and  approximations.  The  whole  procedure  can 
continue  recursively  until  the  output  of  the  filter  reaches  the  minimum  number  of  samples 
equal  to  1.  Thus,  the  maximum  number  of  decomposition  stages  is  log2(N)  for  a  signal  of 
length  N.  Figure  13  shows  a  three-stage  decomposition  tree.  The  original  signal  can  now 
be  written  as  a  combination  of  details  and  approximations  of  various  levels.  A  few 
examples  are: 


•  X=Ai+Di 

•  X=  A3+D3+I)2+Di. 


Note  that  the  detail  Di  occupies  a  large  band  of  the  overall  frequency  band 
when  the  signal  is  decomposed  using  a  two-stage  decomposition  X=A2+D2+Di,  as  shown 
in  Figure  14.  Further,  note  that  the  region  covered  by  Di  is  not  decomposed  again  at  later 
stages  of  the  decomposition,  which  prevents  from  zooming  on  smaller  frequency  bands  in 
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the  region  \fs/4,  f/2]  .  Further  decompositions  on  the  HP  region  are  considered  in  the 
wavelet  packet  analysis,  which  can  be  viewed  as  a  generalization  of  the  wavelet 
transform. 


A 

A2 

D2 

N  / 

V 
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/,\ 

N 
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.•  \ 
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fc/4 

fJ2 

Figure  14:  Spectral  coverage  for  a  two-stage  wavelet  decomposition. 

c)  Wavelet  Packet  Analysis 

The  wavelet  packet  decomposition  is  an  extension  of  the  wavelet 
transform  where  both  details  and  approximations  are  decomposed  further  to  the  higher 
level,  thereby  allowing  flexibility  in  the  decomposition.  This  scheme  leads  to  Nlog2N 
possible  decompositions  for  a  signal  of  length  N,  where  any  complete  level  of  the  tree 
forms  a  complete  orthogonal  basis.  Figure  15  shows  the  decomposition  tree  of  a  wavelet 
packet  at  depth  j=3.  A  few  possible  representations  of  the  signal  are: 

•  X=Ai+AD2+DD2 

•  X=AAA3+DAA3+DA2+AD2+ADD3+DDD3. 
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Figure  15 .'Wavelet packet  decomposition  at  level  3. 


4.  Cosine  Packet  Transfoim 

While  the  wavelet  packet  (WPT)  performs  a  multi-resolution  decomposition  by 
partitioning  the  frequency  axis,  the  cosine  packet  transform  (CPT)  performs  a  multi¬ 
resolution  decomposition  by  partitioning  the  time  axis.  As  a  result,  the  CPT  performs 
better  for  narrowband  signals.  Next,  we  briefly  discuss  the  discrete  cosine  transform 
(DCT)  and  the  local  cosine  transform  (LCT)  and  show  how  they  are  used  to  perform  the 
CPT. 

a)  Discrete  Cosine  Transform  IV 

The  DCT  is  the  inner  product  of  the  signal  with  a  windowed  cosine  (also 
called  blocked  cosine).  Specifically,  the  Discrete  Cosine  Transform-iv  (DCT-IV)  of  a 
signal  x(n)  with  length  N  is  defined  as: 

^  x(n)  •  cos[-^^^  integer.  (01-23) 

I  „  iV 
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Figure  16  shows  two  blocked  cosine  functions  for  N=512  and  k=l  and  2. 


k=1 


Figure  16:  Blocked  cosines,  k=l  and  k=2. 

The  DCT-iv  is  closely  related  to  the  discrete  Fourier  transform  (DFT)  by 
the  following  relation: 

=  (in-24) 

Such  a  relationship  shows  that  the  DCT  is  a  very  fast  algorithm,  as  it  takes 
advantage  of  fast  DFT  algorithms. 

6)  Local  Cosine  Transform 

As  described  earlier,  the  basis  of  the  DCT  is  the  blocked  cosine.  Note  that 
the  rectangular  window  used  in  the  DCT  may  give  undesirable  sidelobes  in  the  frequency 
domain.  A  solution  to  this  problem  is  to  select  a  Gaussian  shaped  window  to  reduce  this 
effect.  In  practice,  the  window  used  has  the  general  form  [7] : 
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r{t)  =  exp(7  •  pit))  ■  sin(^(r)) . 


(in-25) 


The  effect  of  the  window,  also  called  bell,  is  to  localize  the  blocked  cosine 
both  in  time  and  in  frequency  better.  Figure  17  shows  the  local  cosine  for  p(t)=0  and 
(p(t)=7j/4[l  +sin(7Tt)}. 

c)  Cosine  packet  transform 

The  cosine  packet  transform  partitions  the  time  axis  following  the  same 
concept  as  that  used  for  the  frequency  axis  in  the  wavelet  packet  decomposition.  The 
local  cosine  transform  is  applied  to  each  time  segment  of  the  signal. 


Blocked  Cosine 


Figure  17:  Local  cosine.  From  Ref  [4] 

The  resulting  tree  structure  is  the  same  as  that  obtained  with  the  wavelet 
packet.  However,  the  CPT  time  resolution  improves  while  the  frequency  resolution 
worsens,  as  the  decomposition  level  increases. 
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5.  Best  Basis  and  Matching  Pursuit 

Recall  that  we  can  decompose  a  discrete  signal  of  length  N  using  wavelets  or  local 
cosine  packets  with  the  maximum  of  j=log2N  stages.  The  time  frequency  atoms  that 
participate  in  the  full  decomposition  form  an  orthogonal  basis.  Wickerhauser  and 
Coifman  developed  an  algorithm  named  ‘Best  Basis’  which  finds  a  complete  set  of 
orthogonal  basis  functions  that  minimizes  a  user-specified  information  cost  criterion  [7]. 

Another  type  of  decomposition  called  the  matching  pursuit  was  proposed  by 
Mallat  and  Zhang  [8].  The  matching  pursuit  offers  much  more  flexibility  in  the  type  of 
decomposition  than  the  WP  decomposition  does,  as  it  is  not  restricted  to  orthonormal 
basis  functions.  In  this  algorithm,  the  signal  is  decomposed  adaptively  into  a  set  of 
functions,  not  necessarily  orthonormal.  Specifically,  for  given  dictionary  D  =  {gy}y^Y^ 

of  P  vectors  (atoms),  the  algorithm  starts  by  projecting  the  signal  s(n)  onto  each  function 
of  the  dictionary  D  and  computing  the  residue  RS .  Thus, 

S  =  {s,gy^)gy,+RS.  (in-26) 

At  each  iteration,  the  algorithm  selects  the  vector  gy^  with  maximum  inner  product  with 
the  signal  S  and  smallest  ||i?S  ||.  At  the  next  iteration,  the  residue  RS  gets  projected  into 
the  other  functions  of  the  dictionary  to  determine  the  second  decomposition  vector.  The 
procedure  continues  until  the  residue  is  zero  or  its  norm  is  small  enough.  The  main 
drawback  in  this  scheme  is  that  the  decomposition  usually  has  an  error  due  to  the  last 
residue,  unless  the  signal  can  be  decomposed  exactly  using  the  given  dictionary  in  a  finite 
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number  of  steps.  In  addition,  the  computational  load  is  higher  than  that  of  wavelet  packet- 
based  decompositions. 


C.  ENERGY  DISTRIBUTIONS 

So  far  we  discussed  atomic  decompositions  which  lead  to  linear  time-frequency 
representations.  Energy  distributions  distribute  the  energy  of  the  signal  in  both  time  and 
frequency  and  are  viewed  as  quadratic  transformations  of  the  signal.  The  basic  idea  is  that 
the  signal  preserves  its  energy  in  both  time  or  frequency  representation.  According  to 
Parseval’s  relation: 

oo  OO 

=  J  -dt  =  J|S(/)p  -df .  (in-27) 

— oo  — oo 

If  we  consider  the  values  |j(0|^  and  |S(/)|^  as  energy  densities,  then  it  would  be 
convenient  to  find  a  joint  time  frequency  energy  density  Ps(t,f)  such  that  [9]: 

oo  oo 

^5  =  J  fPs(tJ)-dt-df.  (m-28) 

— oo  — oo 

This  chapter  considers  one  such  class  of  distributions:  the  Cohen’s  class. 
Distributions  in  this  class  are  covariant  by  translation  in  time  and  frequency  [10].  The 
general  expression  describing  the  Cohen’s  class  of  energy  distributions  is  given  by: 

oo 

— oo 

(m-29) 
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where  is  called  the  parameterization  function. 


1.  The  Wigner-Ville  Distribution 

The  most  popular  distribution  in  the  Cohen  class  is  the  Wigner-Ville  distribution, 
which  is  defined  as: 

oo 

W,(t,f)  =  +  (m-30) 


or  equivalently  as: 


(ni-31) 


Note  that  equation  (111-30)  is  derived  from  equation  (111-29)  when  =  1 .  The 

Wigner-Ville  distribution  has  many  useful  properties.  For  example,  it  allows  for  perfect 
localization  of  linear  chirps,  with  frequency  variation  defined  asf(t)=at+fo  and  described 


as: 


s{t)  =  e  2 

The  WV  distribution  of  a  linear  chirp  is  given  by: 


(m-32) 


a 


W,(t,f)  =  S(f-(fo+-t)). 


(in-33) 


Further  details  on  the  Wigner-Ville  distribution  properties  may  be  found  in 
references  [2,9].  Note  that  the  main  disadvantage  of  the  WVD  is  the  interference  terms 
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that  appear  when  the  signal  is  not  linearly  modulated  or  more  than  one  signal  is  present. 
These  interference  terms  are  due  to  the  bilinear  nature  of  the  WVD.  Figure  18  shows  the 
WVD  and  the  effect  of  the  cross  terms  on  four  different  signals.  The  effect  is  negligible  in 
the  noise-free  linear  chirp  scenario  only.  Cross  terms  show  as  a  line  between  the  2  true 
components  when  the  signal  is  composed  of  two  parallel  chirps. 


Linear  Chirp  Two  linear  Chirps 


“  Hyperbolic  Chirp 


Figure  18:  Wigner-Ville  distribution,time  is  normalized  over  the  pulse  duration. 

Another  potential  drawback  of  the  WV  distribution  is  the  presence  of  spectral 
aliasing  when  the  signal  is  real  and  sampled  near  its  Nyquist  frequency.  However,  this 
problem  can  be  avoided  by  using  the  analytical  version  of  the  signal  or  by  oversampling  it 
by  at  least  a  factor  of  two. 
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2.  Pseudo  Wigner-Ville  Distribution 

Equation  (111-31)  implies  that  we  must  integrate  from  r  =  -oo  to  which  is 

a  problem  in  practice  [9].  To  overcome  this  problem,  a  windowed  version  of  the  WVD 
called  the  Pseudo  Wigner-Ville  Distribution  (PWVD)  is  defined  as: 

oo 

PW,{t,f)=  +  (in-34) 

— oo 

or  equivalently: 

oo 

PWs(t,f)=  (ni-35) 

— oo 

where  h(t)  is  a  time  window.  Note  that  the  interference  terms  have  a  strong 
oscillatory  component  [2].  Windowing  in  the  time  domain  actually  results  in  some 
frequency  domain  smoothing,  thereby  reducing  the  interferences,  as  shown  in  Figure  19. 
However,  note  that  applying  a  time  window  to  minimize  the  interference  terms  results  in 
worsening  of  the  frequency  resolution.  In  addition,  the  PWVD  looses  many  of  the 
valuable  WVD  properties  [2]. 
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Two  linear  Chirps 


Linear  Chirp  in  Noise  (SNR=2db)  Hyperbolic  Chirp 


Figure  19:  Pseudo  Wigner-Ville  Distribution,  time  is  normalized  over  the  pulse 
duration. 


3.  Smoothed  Pseudo  Wigner-Ville  Distribution 

The  PWVD  provides  an  attenuation  of  the  interference  terms  by  smoothing  in  the 
frequency  domain.  Interference  terms  can  be  further  reduced  by  smoothing  in  the  time 
domain.  The  resulting  distribution  is  called  the  Smoothed  Pseudo  Wigner-Ville 
Distribution  (SPWVD)  and  is  described  as  follows: 

oo  oo 

SPW^{t,f)=  ^h{t)^g{s-T)-s{t  +  Tl2)-s*{t-rl2)-e~^‘^^^dT ,  (III-36) 

—oo  — oo 

where  g(t)  is  a  window  in  the  time  domain.  Note  that  the  SPWVD  becomes  the 
PWVD  when  g(t)=S(t).  The  drawback  of  this  method  is  that  the  frequency  resolution 
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further  degrades.  The  effects  of  the  SPWVD  are  shown  in  Figure  20  which  shows  that  the 
cross  terms  have  almost  been  eliminated  while  the  frequency  resolution  has  been 
degraded. 


Linear  Chirp 


Two  linear  Chirps 


Hyperbolic  Chirp 


Figure  20:  Smoothed  Pseudo  Wigner-Ville  Distribution, time  is  normalized  over 
the  pulse  duration 


4.  The  Reassignment  Method 

We  see  that  the  WVD  resolution  worsens  when  we  try  to  suppress  the  interference 
terms.  The  reassignment  method  was  introduced  in  an  attempt  to  improve  the 
interpretation  of  the  transformation  [9].  This  scheme  assumes  that  the  energy  distribution 
in  the  time-frequency  plane  resembles  a  mass  distribution  and  moves  each  value  of  the 
time-frequency  plane  located  at  a  point  (t,f)  to  another  point  (t’,f ),  which  is  the  center  of 
gravity  of  the  energy  distribution  in  the  area  of  (t,f).  The  result  is  a  very  focused 
representation  with  high  intensity  since  the  value  at  the  point  (t’,f )  is  the  sum  of  all  the 
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neighboring  values.  The  reassignment  method  can  be  applied  to  most  energy 
distributions  as  well  as  to  the  spectrogram.  However,  the  computational  load  is  quite 
high.  Figure  21  shows  the  hyperbolic  chirp  used  in  Figure  18  to  20  as  represented  for  the 
reassigned  spectrogram  distribution,  the  reassigned  PWVD  and  the  reassigned  SPWVD. 
Note  the  quality  of  the  focusing  in  all  three  methods  and  the  lack  of  cross  terms  in  the 
PWVD  and  SPWVD. 


Reassigned  Spectogram 


Figure  21:  Reassignment  method  applied  to  a  hyperbolic  chirp  test  signal,  time  is 
normalized  over  the  pulse  duration. 


IV.  THE  RADON  TRANSFORM 


Earlier  chapter  presented  several  techniques  to  construct  an  image  representation 
of  a  signal  in  the  time-frequency  plane.  The  next  step  is  to  extract  the  information  from 
this  image  to  recover  the  modulation  parameters  for  both  linear  and  hyperbolic  chirps.  An 
important  tool  in  image  processing  is  the  Radon  transform  which  is  used  to  extract  line 
information.  We  will  briefly  review  and  then  apply  this  transform  to  extract  linear  chirp 
parameters  [13]. 

A.  INTRODUCTION 

Let’s  assume  we  have  an  arbitrary  function  f(x,y)  in  a  subspace  of  R  .  The  two- 
dimensional  Radon  transform  is  defined  as  the  projection  or  line  integral  of  the  function 
f(x,y)  along  all  possible  lines  L  [13].  Mathematically,  the  transformation  is  described  as: 

R  =  J  f(x(s,L),  y(s,L))ds .  (IV-1) 

L 

Recall  that  the  equation  of  a  line  in  polar  coordinates  is  given  by: 

p  =  X- cos(t9) -I-  y  •  sin(t9) ,  (IV -2) 

where  p  and  t?  represent  the  distance  from  the  origin  and  the  angle  measured 
counterclockwise  from  the  x  axis  respectively,  as  shown  in  Figure  22.  Now  suppose  we 
use  another  coordinate  system  with  axes  rotated  by  the  angle  r?.  The  new  x-axis  lies  on 
the  line  with  associated  orthogonal  direction  “s”.  The  two  cartesian  systems  xoy  and  pos 
are  related  to  each  other  via  the  following  relation: 
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(IV-3) 


cost? 

sint? 

X 

s 

-sint? 

cost? 

'  y_ 

Equation  (IV-1)  can  be  rewritten  as: 


-  5  sin  t?,  yO  sin  t?  +  j  COS  t?)  • 


(IV-4) 


The  above  equation  shows  that  the  Radon  transform  translates  a  two-dimensional 
function  of  the  variables  (x,y)  to  one  with  variables  .  Thus,  the  Radon  transform  of 
an  image  taken  at  a  specific  angle  &  is  the  projection  of  the  image  onto  the  line  which 
forms  an  angle  ^  with  the  x-axis. 


Figure  22:  Two-  dimensional  Radon  transform. 

The  Radon  transform  of  a  single  line  with  a  slope  angle  ^  for  the  specific  angle 
^  =  Y  +  ^  is  a  single  point  with  intensity  equal  to  the  sum  of  the  intensities  at  each  point 
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on  the  line.  This  property  allows  detection  of  lines  in  an  image.  In  addition,  the  transform 
is  also  robust  to  noise  degradations. 

B.  LINE  PARAMETER  roENTIFIGATIONS 

Line  parameters  can  be  obtained  using  the  Radon  transform  as  follows:  Assume 

the  image  under  investigation  contains  a  line  which  has  an  angle  (p  with  the  x-axis,  as 
shown  in  Figure  23.  The  equation  of  the  line  can  be  defined  in  terms  of  its  slope  a  and 
initial  offset  value  b. 

y  =  a-  x-\-b,  a  =  tan(^) .  (TV-S) 

The  Radon  transform  RCp.d)  of  the  image  for  angles  0°  to  180°  is  maximum  when 
the  projection  of  the  line  has  a  minimum  area.  Thus,  it  is  maximum  when  the  line  is 
perpendicular  to  the  projection  line,  i.e.,  when  ‘d=(p+9(f,  which  leads  to 

a  =  tan(t?-90°).  (IV-6) 

Now,  if  we  take  the  Radon  transform  of  the  image  for  the  specific  angle 
■&=(p^9(f,  we  can  estimate  the  offset  parameter  b  from  the  position  of  its  maximum  value 
along  the  axis  p.  This  position  is  indicated  in  Figure  23  as  “C”.  Next,  the  offset  parameter 
b  can  be  computed  as: 

OC— ^sin(t9-90°)  ^  OC  +  ^cos(j9) 

=  + - 2 - =  il2L  + - 2 - ^  (IV.7) 

2  cos(?9-90°)  2  sin(r?) 

where  Ny  is  the  length  of  the  vertical  side  of  the  image,  while  N  is  the  length  of 
the  horizontal  side  (in  this  case  Nx=Ny=256).  Note  that  the  distance  OC  can  be  positive  or 
negative  and  is  actually  negative  in  Figure  23. 
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Thus,  the  equation  of  the  line  can  be  determined  by  first  applying  the  Radon 
transform  of  the  image  for  angles  between  0°  to  180°,  then  finding  the  coordinates  of  the 
maximum  point  Pmax),  as  shown  in  Figure  24,  and  finally  using  equations  (IV-6) 

and  (IV-7)  using  OC=Pmax  and  Note  that  the  accuracy  of  the  estimates  depends 

on  the  resolution  of  the  image  and  the  size  of  the  angle  step  selected  for  the  Radon 
transform.  Unfortunately,  any  attempt  to  increase  the  resolution  (size)  or  decrease  the 
angle  step  results  in  increasing  the  computation  load. 


Figure  23:  Geometry  of  the  Radon  transform  for  an  image  containing  a  single 
line. 


50  100  150 

Angle  in  degrees 


Figure  24:  Radon  transform  for  the  line  shown  in  Figure  23  for  &=(f...l80°,  f 
increment. 
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To  minimize  the  computational  load  we  apply  the  Radon  transform  in  two  stages. 
First,  we  scan  the  image  from  0°  to  179°  in  steps  of  2°.  We  determine  the  angle  dim  that 
corresponds  to  the  column  that  includes  the  point  with  the  highest  intensity.  Next,  we 
scan  the  image  from  the  angle  dim-l°  to  dim+l°  in  steps  of  0.1°.  Using  this 
implementation  allows  the  maximum  line  slope  error  to  decrease  to  ±0.05°  without 
having  to  cover  the  whole  range  of  angles  at  that  fate. 
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V.  LINEARLY  MODULATED  CfflRPS 


This  chapter  presents  the  scheme  used  to  estimate  linear  chirp  parameters.  This 
chapter  is  divided  in  four  sections.  The  first  part  presents  the  method  used  to  generate  the 
signals.  The  second  section  presents  the  time-frequency  transformations  considered  in 
this  study  from  which  the  linear  chirp  parameters  are  extracted.  The  third  section 
describes  the  image  processing  technique  applied  to  extract  the  features  from  the  time 
frequency  image.  Finally,  the  last  section  compares  the  results  obtained  for  each  time 
frequency  representation  considered  and  their  robusmess  to  noise  degradations. 

A.  SIGNAL  GENERATION 

The  radar  signals  considered  in  this  study  are  synthetic  and  consist  of  a  train  of 
several  linearly  modulated  pulses,  as  shown  in  Figure  25. 


t 


Figure  25:  Synthetic  radar  signal,  linear  chirp  modulation. 
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First,  we  assume  that  we  can  isolate  one  single  pulse  of  duration  x.  This  extraction 
can  easily  be  done  with  an  energy  detector  in  medium  to  high  SNR  levels.  Thus,  the 
received  signal  has  an  instantaneous  frequency  f(t)  defined  as: 

/(0  =  /o+^-?  (V-1) 

Some  of  the  time-frequency  transformations  considered  in  the  study  use  the 
analytical  version  derived  from  the  real  received  signal.  In  such  cases,  the  analytical 
signal  is  computed  with  the  Hilbert  transform.  The  analytical  noise-free  linear  chirp 
signal  s(t)  is  given  as: 


s(t)  =  e 


(V-2) 


The  signal  is  assumed  to  be  sampled  a  rate  of  512  [samples/pulse  length].  Recall  that  the 
sampling  frequency  for  a  real  signal  must  be  equal  to  at  least  twice  its  highest  frequency 
component.  Thus,  the  signal  may  be  heterodyned  down  to  a  lower  frequency  if  needed. 
The  discrete  real  part  of  the  signal  is  given  as: 

s[n]  =  Re(e  •  )  =  cos(2;r(^-n  +  ^.n^)),  n=l,..,512  (V-3) 

J s  s 


We  can  rewrite  the  above  equation  as: 


s[n]  =  cosCln{f^Q-n  +  ^-n^))  where  fno=^,  a  =-^  , n=l,..,512 

2  fs 


(V-4) 
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The  terms  a  and  /,o  represent  the  normalized  slope  and  the  normalized  starting 
frequency  respectively,  with  respect  to  the  sampling  frequency  fs.  Thus,  the  normalized 
frequency  equation  for  the  linear  chirp  discrete  signal  is  given  by: 

/„=/„o+a-n  ,  n=l,..,512.  (V-5) 

Multiple  linear  chirp  signals  trials  were  generated  by  randomly  selecting  both  the 
initial  frequency  and  the  slope.  Note  that  sampling  constraints  need  to  be  satisfied  to 
avoid  aliasing  in  the  resulting  discrete  linear  chirp.  For  example,  the  parameter  a  needs  to 
be  selected  so  that  the  final  instantaneous  frequency  doesn’t  exceed  0.5  for  a  given  initial 
frequency  ^0.  Thus,  aliasing  may  be  avoided  by  selecting  both  the  initial  and  final 
normalized  instantaneous  frequencies  within  the  range  [0,  0.5],  and  computing  the 
corresponding  slope  parameter  a.  Initial  and  final  frequencies  are  selected  from  a 
uniform  distribution  defined  in  the  range  [0,0.5].  Next,  the  chirp  is  corrupted  by 
additive  white  Gaussian  noise.  Analytical  expressions  needed  to  compute  some  of  the 
time-frequency  transformations  are  obtained  with  the  Hilbert  transform  of  the  noisy 
signal.  Real  signal  expressions  were  used  to  compute  the  wavelet-based  decompositions. 

B.  SIMULATION  SET  UP  AND  FEATURE  EXTRACTION 

Eleven  different  time-frequency  and  wavelet-based  representations  were 

considered  in  this  study.  The  goal  was  to  select  a  small  number  of  transformations 

leading  to  the  best  “image  quality”  from  that  set.  The  representations  considered  were: 

-Wavelet  packet  best  basis 
-Cosine  packet  best  basis 
-Wavelet  pursuit 

43 


-Cosine  pursuit 
-Wigner-Ville  distribution 
-Pseudo  Wigner-Ville  distribution 
-Reassigned  Pseudo  Wigner-Ville  distribution 
-Smoothed  Pseudo  Wigner-Ville  distribution 
-Reassigned  Smoothed  Pseudo  Wigner-Ville  distribution 
-Spectogram  distribution 
-Reassigned  spectogram  distribution 

Wavelet  and  cosine  decompositions  were  implemented  with  the  software 
“Wavelab  7.01”  software  [12],  while  all  others  used  the  ‘Time-Frequency  Toolbox” 
[11].  Figure  26  shows  the  resulting  images  obtained  for  a  linear  chirp  with  SNR=10dB, 
when  no  energy  thresholding  is  applied  to  the  images.  A  few  comments  are  in  order 
regarding  the  selection  of  the  transform  parameters. 

1.  Number  of  Atoms 

Recall  that  defining  a  line  in  a  plane  requires  two  points  only.  However,  these  two 
points  must  be  perfectly  localized  in  both  time  and  frequency,  and  must  be  inunune  from 
any  noise  interference.  Theoretically,  we  could  use  two  atoms  only  from  an  atomic 
decomposition  to  define  the  linear  frequency  equation.  However,  this  doesn’t  hold  in 
practice  as  the  atoms  are  not  perfectly  localized.  Note  that  a  larger  number  of  atoms  may 
better  represent  the  line  trend.  However,  some  of  the  atoms  may  represent  noise 
contributions  for  noisy  signals.  Therefore,  selecting  the  number  of  atoms  to  represent  a 
linear  chirp  in  noise  requires  a  trade-off  between  these  two  issues:  fewer  atoms  to  denoise 
the  signal  and  more  atoms  to  improve  the  resolution.  We  selected  10  atoms  per 
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decomposition  for  the  atomic  decompositions  used  in  our  study  after  running  several  trial 
cases. 

2.  Maximum  Decomposition  Level 

Next,  the  maximum  decomposition  level  was  set  to  7,  as  simulations  showed  no 
advantage  in  going  to  higher  levels. 

3.  Wavelet  Type 

The  mother  wavelet  function  was  selected  after  several  trials  among  the  readily 
provided  functions  in  the  Wavelab  software  [12].  Finally,  we  selected  Daubechie-IV 
since  it  gave  the  best  time-frequency  representation  for  the  types  of  signal  considered. 

4.  Image  Thresholding 

No  image  intensity  thresholding  was  applied  to  the  time-frequency  images,  as 
simulations  showed  that  it  this  step  worsened  the  results  for  low  SNR  levels.  Note  that 
implicit  denoising  is  actually  performed  by  selecting  a  small  number  of  atoms  for  the 
atomic  decompositions. 

5.  Window  Length 

The  PWV,  SPWV  and  Reassigned  SPWVD  transformations  use  a  frequency 
window  which  has  a  Hamming  (N/4)  time  domain  expression,  and  a  Hamming  (N/10)  for 
time  smoothing  window.  The  analysis  window  selected  for  the  spectrogram  is  Hamming 
(N/4),  where  N  is  the  length  of  the  signal  (512  bins). 
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Wavelet  Packet 


Cosine  Packet 


Figure  26:  Various  time  frequency  representations,  linear  chirp,  SNR=10dB. 

6.  Radon  Transformation  Implementation  Issues 

The  Radon  transform  is  selected  to  extract  the  line  equation  from  the  time- 
frequency  image,  as  described  earlier  in  Section  lH.  We  use  a  two-stage  implementation 
with  final  degree  increment  0.1°.  The  size  of  the  resulting  image  for  each  time  frequency 
representation  is  set  to  256x256  points,  as  the  Radon  transform  for  larger  images  would 
be  too  computationally  expensive  for  a  MATLAB  implementation.  One  hundred 
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randomly  generated  linear  chirps  for  a  given  SNR  level  were  generated,  and  the  SNR 
varied  in  the  range  -lOdB  to  lOdB.  Next,  the  chirp  parameters  were  estimated  from  the 
images  generated  by  each  of  the  eleven  time-frequency  transformations  considered. 
Performance  comparisons  are  presented  next. 

C.  PERFORMANCE  COMPARISONS 

1.  Evaluation  of  Time-frequency  Representations 
Recall  that  the  maximum  theoretical  slope  error  is  ±0.05°  since  the  final  step 
angle  in  the  radon  transform  is  0.1°,  as  discussed  in  Section  IV.  However,  we  also  have  to 
add  noise  effects,  and  quantization  errors  introduced  by  the  finite  resolution  in  the  image. 
Figures  27  to  29  present  the  mean  and  the  standard  deviation  of  the  absolute  slope  and 
offset  errors  as  a  function  of  the  SNR  level  for  all  time-frequency  transformations 
considered  in  the  study. 
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Figure  28:  Linear  chirp  signal;  Slope  and  offset  error  for  the  PWV,  reassigned 
PWV,  SPWV  and  reassigned  SPWV  decompositions  (SNR  in  dB.) 


Degrees 


SNR  SNR 

Figure  29:  Linear  chirp  signal;  Slope  and  offset  error  for  the  Wigner-Ville 

decomposition,  Spectrogram,  and  reassigned  Spectrogram  (SNR  in  dB.) 

A  few  comments  are  in  order: 

1.  Results  show  that  the  energy  distributions  perform  better  than  the  atomic 
decompositions.  This  is  to  be  expected  as  they  provide  a  more  accurate 
“image”  in  the  time-frequency  plane.  Note  that  atoms  cannot  be  well 
localized  in  both  time  and  frequency,  as  would  be  required  to  represent 
linear  chirps  accurately.  The  best-basis  cosine  packet  decomposition  gives  the 
best  results  followed  by  the  cosine  pursuit  scheme  for  atomic  decompositions. 
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2.  Most  of  the  energy  distributions  have  slope  errors  very  close  to  the  theoretical 
value  of  0.05°  for  medium  to  high  SNR  levels. 

3.  The  Wigner-Ville  distribution  has  a  very  good  and  almost  stable  performance 
for  SNR’s  in  the  range  of  -6dB  to  lOdB.  The  smoothing  time  window  present 
in  the  Pseudo  Wigner-Ville  distribution  improves  the  estimation  in  higher 
SNR  but  shrinks  the  effective  range.  The  Smoothed  Pseudo  Wigner-Ville 
distribution  has  a  slightly  worst  performance  at  high  SNR  but  also  has  the 
widest  effective  range.  The  smoothing  in  time  and  in  frequency  eliminates  the 
interference  terms  almost  completely  so  that  the  representation  is  more 
immune  to  the  noise  than  other  transformations  are.  However,  frequency 
smoothing  results  in  lowered  performance  at  high  SNR  levels. 

4.  Results  show  that  the  reassign  method  usually  improves  the  performance  at 
high  SNR  levels,  as  it  forces  the  representations  to  be  more  “focused”. 
Unfortunately,  applying  the  reassign  method  in  low  SNR  levels  worsens  the 
performances.  This  is  to  be  expected  as  the  presence  of  noise  close  to  the  line 
moves  the  local  center  of  gravity  of  the  distribution  away  from  its  theoretical 
value. 

Table  1  lists  the  mean  absolute  error  for  the  frequency  slope  and  offset  parameters 
at  SNR  equal  to  lOdB  for  all  transformations  considered  for  a  “rough”  comparison  of  the 
transformation  performances.  In  addition,  Table  1  lists  the  SNR  level  at  which  the  errors 
suddenly  increase,  indicating  the  SNR  level  at  which  using  a  specific  distribution  may 
become  more  questionable. 
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Mean  of  absolute  slope 
error  @  SNR=10dB 
(in  degrees) 

Mean  of  absolute  offset 
error  @  SNR=10dB 
(in  bins) 

Breaking  point 
(in  dB) 

Wavelet  Packet 

0.375 

4 

0 

Cosine  Packet 

0.5149 

3.98 

0 

Wavelet  Pursuit 

0.37 

2.35 

4 

Cosine.  Pursuit 

0.4297 

2.35 

0 

WV 

0.1198 

1.08 

-7 

PWV 

0.0623 

1.02 

-4 

Reassigned  PWV 

0.0729 

1 

-5 

SPWV 

0.0792 

1.02 

-7 

Reassigned 

SPWV 

0.0563 

1.02 

-7 

Spectrogram 

0.098 

1.77 

-7 

Reassigned 

Spectrogram 

0.12 

1.78 

-7 

Table  1:  Mean  absolute  errors  for  frequency  slope  and  offset.  100  trials, 
SNR=10dB,  performance  breaking  points  for  all  time-frequency  methods 
considered. 


2.  Error  Evaluation 

Recall  the  received  signal  was  sampled  at  a  rate  of  512  [samples/pulse  period] 
while  the  image  size  was  set  to  256x256bins  in  order  to  speed  up  the  computations. 
However,  the  errors  presented  were  corrected  to  correspond  to  an  image  with  size 
512x256,  where  the  time  and  frequency  axes  have  512  and  256  bins  respectively.  Next, 
the  goal  is  to  investigate  how  the  estimated  errors  relate  to  the  frequency  error  obtained 
for  the  frequency  expression  given  in  equation  (V-1). 
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Recall  that  the  normalized  frequency  expression  is  given  by: 

/« 

and  the  slope  angle  d  measured  in  the  image  is: 

a-N-2Nf 

tan(z?)  = - — — (V-6) 

where  a  is  the  frequency  slope  parameter,  N  is  the  length  of  the  signal,  and  Nf  and 
Nt  represent  the  number  of  bins  in  the  frequency  and  time  axis  respectively.  In  our  case 
Nf=N/2  and  Nt=N.  Thus,  equation  (V-6)  simplifies  to: 

tan(i?)  —  a- N .  (V-7) 

Recall  that 

a  f  k 

s[n]  =  cosi2n:(f„o-n  +  --n^)),  where  a  =— , n=l,..,512 

^  f  s  f s 

Let’s  assume  that  the  linear  chirp  has  frequency  slope  ko  and  that  is  the 
corresponding  angle,  which  leads  to: 

=  .  (V-8) 

fs  ^ 

Now,  let’s  assume  we  estimate  an  angle  equal  to^^j/  instead  of  ^  due  to  errors  in 
the  estimation  process.  Thus, 

(V-9) 

where  Afr  is  the  angle  error.  As  a  result,  the  resulting  estimated  frequency  slope 
kest  expression  is  given  by: 
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Equations  (V-11)  and  (V-15)  provide  the  relations  between  errors  measured  in  the  time- 
frequency  image  representation  of  size  N/2xN,  where  N  is  the  length  of  the  signal,  and 
the  corresponding  errors  in  the  initial  time-varying  frequency  expression  given  in 
Equation  (V-1).  Note  that  Afb  is  measured  in  bins. 

3.  Multi-pulse  Processing 

Up  to  this  point,  all  results  were  derived  by  extracting  the  frequency  modulation 
parameters  from  one  noisy  radar  pulse.  Performances  improve  when  using  more  than  one 
pulse.  Assume  we  isolate  five  radar  pulses,  where  each  pulse  contains  the  same 
transmitted  signal  at  a  given  SNR  level.  Note  that  the  signal  information  gets  mapped  to 
the  same  location  of  the  time-frequency  plane,  while  the  noise  contributions  scatter  to 
different  location  in  each  trial.  Thus,  averaging  multiple  time-frequency  images  improves 
the  quality  of  the  signal  information. 

Thus,  averaging  the  contribution  of  five  pulses  leads  to  the  results  shown  in 
Figure  30  obtained  for  the  Smoothed  Pseudo  Wigner-Ville  transformation.  One  hundred 
realizations  per  SNR  level  are  used,  and  SNR  levels  varied  from  -20  to  OdB.  Note  that 
there  is  no  need  to  consider  higher  SNR  levels  as  the  error  curves  already  flatten  for 
SNR=0dB.  Results  show  a  significant  improvement  over  using  one  pulse  only.  Further 
improvements  may  be  obtained  by  considering  a  larger  number  of  pulses.  However,  this 
will  result  in  a  direct  increase  of  the  computation  time. 
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Figure  30:  Slope  and  offset  error  of  SPWV  distribution;  five  realizations. 

4.  Examples 

We  apply  the  above  results  to  a  specific  radar  using  pulse  compression  readily 
available  to  us:  the  SPS-40B  radar.  The  Radar  Set  AN/SPS-40B  is  used  for  the  detection 
ranging  and  tracking  of  air  targets  at  long  range  in  American  and  foreign  navy  services 
[24].  In  this  pulse  mode,  the  pulse  length  is  x=60psec.  The  transmitted  waveform  inside 
the  pulse  is  a  linear  modulated  down-chirp  with  a  bandwidth  of  2MHz  centered  at  a 
frequency  that  can  be  varied  manually  by  the  operator  from  402.5  to  447.5  MHz. 

Let’s  assume  that  the  radar  operates  at  437.5  MHz,  i.e.,  the  chirp  starts  with  a 
frequency  438.5MHz  and  decays  linearly  to  436.5  MHz.  Our  ESM  (Electronic  Support 
Measures)  receives  a  series  of  noisy  pulses  from  this  radar.  Further,  assume  we  can 
isolate  one  pulse  perfectly  and  that  the  estimated  received  signal  is  centered  at  438MHz 
with  a  bandwidth  equal  to  2MHz.  Given  that  the  duration  of  the  pulse  is  bOpsec  and  if  we 
use  512  [samples/pulse  duration],  the  sampling  frequency  should  be  set  at: 
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(V-16) 


/.  =  —MHz  =  Z.5?>MHz . 

This  sampling  frequency  does  not  satisfy  the  Nyquist  rate.  Thus,  we  can  either 
increase  the  sampling  frequency  or  heterodyne  to  a  lower  center  frequency.  The  first 
option  has  the  drawback  of  increasing  the  number  of  samples  to  deal  with,  and  the 
computational  load  of  the  estimation  schemes  considered.  Heterodyning  the  signal  down 
to  the  baseband  can  be  accomplished  by  multiplying  the  received  signal  with  a  cosine 
function,  and  using  a  lowpass  filter  to  extract  the  information.  Assume  the  heterodyning 
frequency  is  selected  as  440MHz.  Thus,  the  resulting  chirp  signal  is  an  up-chirp  at  the 
frequency  440-437.5=2.5MHz  with  bandwidth  equal  to  2MHz  after  heterodyning  and 
filtering.  The  theoretical  time-frequency  representation  is  shown  in  Figure  31.  Thus,  the 
instantaneous  frequency  of  the  heterodyned  signal  is  given  by: 

/(0  =  1.5-10® +3.33-10'° -r.  (V-17) 

The  heterodyned  signal  may  be  sampled  with  a  frequency  fs=8.53Mhz, 
corresponding  to  5 1 2  samples/pulse  duration. 

Assume  we  use  the  SPWV  distribution  and  the  chirp  parameter  estimation 
procedure  described  above  in  Section  IV-B.  Figure  28  shows  that  the  mean  slope  error  is 
0.09°  and  the  mean  offset  error  is  1  bin  when  the  SNR  is  equal  to  5dB.  Using  equations 
(V-11)  and  (V-15)  leads  to: 

|Ak|=2.610*  and  |Afo|=1.6710^ , 
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meaning  that  the  normalized  errors  for  the  frequency  slope  and  offset  are 

,  M  |A/nl 

equal  to  —  —  0.008  or  0.8%  and - —  =  0.0104  or  1.04%  respectively 

3.33-10*°  1.2-10^ 

for  the  heterodyned  signal. 

Note  that  the  normalized  slope  error  remains  the  same  for  the  original 
received  signal,  as  the  signal  bandwidth  has  not  been  changed  by  the  heterodyning 
process,  while  the  normalized  frequency  offset  error  becomes  equal  to 
1.6710'‘/438.510^  i.e.,  0.0038%. 


Figure  31:  Instantaneous  frequency  for  the  heterodyned  signal. 


58 


VI.  HYPERBOLIC  MODULATED  CHIRPS 

This  section  considers  the  estimation  of  hyperbolic  chirp  parameters.  As  before, 
the  starting  point  is  the  time-frequency  image  representation  of  the  information. 
However,  generalizations  of  the  Radon  transform  to  hyperbolic  line  types  were  not  used. 
Instead,  we  consider  an  iterative  procedure  to  estimate  the  chirp  parameters. 

Note  that  we  assume  we  know  the  type  of  modulation  transmitted,  as  we  do  not 
address  issues  related  to  distinguishing  between  linear  and  hyperbolic  modulation  types. 
Such  classification  issues  are  left  for  future  work. 

A.  SIGNAL  GENERATION 
1.  Introduction 

Assume  that  we  can  isolate  individual  pulses  of  duration  t=512  samples.  Thus,  a 
received  signal  with  hyperbolic  modulation  frequency  is  given  by: 

jc(0  =  cos(2;r(Aln(t  +  ro)  +  ®’0)j  (YI-1) 

where 

/(r)  =— +fi.  (VI-2) 

t+tQ 

The  analytical  signal  s(t)  obtained  from  x(t)  with  a  Hilbert  transform  is  given  by: 

4r)  =  _  (VI.3) 

The  corresponding  discrete  signal  \[n]  is  given  by: 

;2ff(A-ln(n+n<,)-Aln(/,)+— n) 

x[n\  =  e  fs  f,  fs  ^  (VI-4) 
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or  finally 


where: 


x[n]  =  K‘e 


y‘2-;r(fl-ln(n+«Q  )+Z?n) 


(VI-5) 


a  =  A,  b  =  ^,  nQ=f^-tQ,andK  =  . 

/  S 

The  normalized  instantaneous  frequency  expression  is  given  by: 

/„W=— ^+i>. 

n  +  riQ 


(VI-6) 


(VI-7) 


Recall  that  we  need  to  select  the  parameters  a,  b  and  no  such  that  the  range  of  the 

normalized  frequency /,/n7  is  between  0  and  0.5  to  prevent  aliasing,  which  leads  to  the 
following  ranges: 


ae[0,oo],  «o  ^  [0,oo],  fee  [-00,0.5]  (VI-8) 

A  valid  selection  for  a,  b  and  no  is  obtained  by  selecting  one  parameter  in  its 
allowed  range,  and  then  the  other  two  so  that  they  also  fall  in  their  allowed  ranges.  In  the 
simulation,  we  first  select  a  value  for  b  from  a  uniform  distribution  in  the  region  [0, 0.5]. 
The  starting  frequency  at  n=0,  fn(0)  ,is  selected  from  a  uniform  distribution  in  the  region 
[b,  0.5].  The  final  frequency /„(iy),  for  n=N,  where  N  is  the  signal  length,  can  be  selected 
in  the  range  [b,fn(0)].  Defining  b,  fn(0),  and  f„(N)  leads  to  the  following  values  for  a  and 
no: 


^_ifAO)-b)-{b  +  f„m-N 
(/„(0)-/„(iV)) 


(VI-9) 
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(VI-10) 


^  ib  +  f,(N))-N 

®  (fniO)-fn(N))' 

Note  that  this  random  selection  can  sometimes  lead  to  a  frequency  equation  with  a 
very  steep  slope  at  the  beginning,  as  shown  in  Figure  32.  Such  a  behavior  is  undesirable 
because  none  of  the  time-frequency  representations  gives  good  resolution  in  the  area  of 
the  steep  slope,  and  this  modulation  type  is  not  typical  in  radar  applications. 


Theoretical  Time-Frequency  diagram 


n 


Figure  32:  Normalized  frequency  fora=10,  no=25,  b=0.05. 

2.  Hyperbolic  Line  Parameter  Constraints 

We  wish  to  restrict  our  random  selection  of  the  signal  parameters  to  those  leading 
to  chirps  without  steep  slopes,  and  need  to  automate  the  process  so  that  we  can  run 
multiple  trials  to  study  the  scheme  robustness  to  noise  degradations.  Thus,  we  define  a 
figure  of  merit  that  describes  the  amount  of  curvature  present  in  a  given  hyperbolic  line. 


61 


The  parameter  selected  to  characterize  the  curvature  of  the  chirp  is  the  distance  d  defined 
as  the  maximum  distance  between  any  point  of  the  curve  and  its  cord. 

Note  the  parameter  b  shifts  the  hyperbolic  chirp  up  or  down  without  affecting  the 
shape,  thereby  its  curvature.  Therefore,  we  assume  that  b=0  for  simplicity.  In  such  a 

case,  the  chord  associated  to  the  hyperbolic  chirp  /„[«]  =  — - — is  a  line  passing  through 

n  +  n© 


the  points  (0,—)  and  (N , — - — )  with  equation: 

Uq  N  +  nQ 

y.__^  a  _ 

_ np  _  N  +  Hq  Wq 

n-0  N-0 

The  resulting  chord  equation  is  given  by: 

a-n  +  riQ-iN  +  nQ)- f  -a-{N  +  nQ)  =  0,  (VI-11) 

which  leads  to: 


/  = 


-a  a 

- Tl  H - . 

np-CV  +  np)  np 


(VI-12) 


Recall  that  the  gradient  of  a  curve  at  a  given  point  is  a  line  that  passes  through  that 
point  and  has  a  slope  equal  to  the  derivative  of  the  curve  at  that  point.  Thus,  the  gradient 
at  any  point  of  the  hyperbolic  line  is: 


A  =  /'  = 


-a 

(n-t-rap)^ 


(VI-13) 
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Figure  33:  Maximum  distance  between  a  hyperbolic  line  and  its  chord.  a=12, 
n0=25,  b-0. 


The  distance  d  is  maximum  at  the  point  where  its  gradient  is  parallel  to  the 

chord.  As  a  result,  the  gradient  slope  at  the  point  n=^is  equal  to  the  slope  of  the  chord. 

Thus  the  position  of  the  point  n=|can  be  estimated  from  (VI-11)  and  (V-13)  as: 

-a  _  -a 
ii  +  nof  «o‘(^  +  "o)’ 


which  yields  the  coordinates  of  the  point  ^  as: 


^  =  V"o(^  +  ”o)-”0’ 


(VI-14) 


/(!)  = 


a 


^nQ(N  +  nQ) 


(VI-15) 


Recall  the  distance  d  between  a  point  (xi,yi)  and  a  line  with  equation 
•  X  +  ^2  •  y  +  =  0  is  given  by  the  equation: 
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(VI-16) 


j  _\a\-xi+a2-y^+a-i\ 

I  2  2” 

V^i  +^2 

Thus,  the  distance  between  the  chord  and  the  point  can  be  computed  using 

equations  (VI-11),  (VI-14),  (VI-15),  and  (VI-16),  which  leads  to: 

^_|a-^  +  «o-(iV  +  no)-/(|)-a-(V  +  no)| 

+no(V  +  no)^ 

and: 


d  = 


|a|-  (N  +  2-no)-2-^no(-^  +  ”o) 


(VI-17) 


■y/a^+no(A/’  +  no)^ 

The  distance  d  can  be  viewed  as  a  figure  of  merit  for  the  curvature  of  the 
hyperbolic  line.  A  high  value  of  d  means  the  curvature  is  high  and  the  time-frequency 
representation  near  the  time  origin  is  poor.  However,  very  small  values  of  d  represent 
cases  for  which  the  hyperbolic  curve  is  almost  indistinguishable  from  a  straight  line. 
Thus,  we  restricted  the  chirp  signals  generated  so  that  the  distance  d  is  in  the  region 


^  f  80 

'  Js  ■>  J  5 


1024 


1024 


to  avoid  such  cases. 


Next,  additive  zero-mean  white  Gaussian  noise  is  added  to  generate  the  noisy 
chirp  with  a  specific  SNR  level. 


B.  FEATURE  EXTRACTION 
1.  Introduction 

The  signal  time-frequency  representation  can  be  obtained  with  any  of  the  energy 
distributions  discussed  earlier  in  Section  HI.  However,  note  that  the  atomic 
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decompositions  are  not  as  well  suited  as  the  energy  distributions  for  the  task  at  hand,  as 
they  do  not  describe  the  hyperbolic  line  curvature  accurately.  Therefore,  we  only  consider 
energy  distributions  in  this  chapter. 

At  this  point  the  task  is  to  extract  the  three  unknowns  parameters  (a,b,no)  for  a 
given  time-frequency  representation.  The  basic  Radon  transform  can  no  longer  be 
applied,  as  it  is  defined  for  straight  lines  only.  The  Radon  transform  was  extended  to 
detect  functions  of  arbitrary  shape  [17],  however  the  computational  load  is  significantly 
higher. 


Figure  34:  Extraction  of  the  instantaneous  normalized  frequency  expression  from 
the  time-frequency  image.  Before  and  after  median  smoothing,  L=5. 
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The  method  considered  here  approaches  the  problem  quite  differently.  It  is  an 
adaptive  procedure  which  fixes  one  of  the  three  unknown  parameters  at  each  iteration, 
and  estimates  the  other  two.  The  resulting  scheme  is  presented  next. 

2.  Method  Description 

The  instantaneous  frequency  expression  for  the  hyperbolic  chirp  is  extracted  from 
the  two-dimensional  image  by  selecting  the  peak  values  obtained  at  each  time  bin  of  the 
image.  The  result  is  a  vector  containing  the  frequency  values  for  each  time  bin,  as  shown 
in  Figure  34. 

Note  that  the  2-D  instantaneous  frequency  approximation  is  very  accurate  at  high 
SNR  levels,  and  degrades  as  the  SNR  level  decreases,  hi  noisy  environments,  the  pixel 
with  the  highest  energy  obtained  at  a  given  time  bin  may  be  an  outlier,  resulting  in  spikes 
in  the  instantaneous  frequency  estimate,  as  illustrated  in  the  middle  plot  of  Figure  34. 
Such  outliers  can  be  smoothed  out  with  a  median  filter.  Simulations  showed  that  a 
median  filter  of  length  5  smoothed  out  potential  “spikes”  without  loss  of  resolution. 

We  selected  the  three  energy  transformations  leading  to  the  best  time-frequency 
image  quality  for  the  linear  chirp  case,  and  restricted  our  hyperbolic  chirp  analysis  to 
those.  The  distributions  selected  are: 

-  Pseudo  Wigner-Ville  distribution, 

-  Reassigned  Pseudo  Wigner-Ville  distribution, 

-  Reassigned  Smoothed  Pseudo  Wigner-Ville  distribution. 

We  set  the  dimension  of  the  image  at  512x512  bins  to  increase  the  image 
resolution  and  reduce  quantization  errors.  In  addition,  simulations  showed  that  these  three 
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energy  distributions  do  not  necessarily  perform  well  at  the  beginning  and  end  of  the 
image.  As  a  result,  we  only  consider  the  image  from  time  bin  60  to  time  bin  450. 

The  second  step  in  the  proposed  scheme  approximates  the  instantaneous 
frequency  adaptively  with  a  hyperbolic  line  by  minimizing  the  squared  error  between  the 
information  contained  in  the  image  and  a  theoretical  hyperbolic  curve  expression.  Note 
that  the  problem  to  be  solved  is  non-linear,  due  to  the  specific  frequency  law  to  be 
approximated,  as  we  estimate  the  parameters  a,  b  and  no  given  in  equation  (VI-7).  We 
first  tried  to  solve  the  problem  using  a  classical  nonlinear  least  square  iteration  scheme 
provided  with  the  function  “Isqnonlin”  from  the  MATLAB  optimization  toolbox  [15]. 
Simulations  showed  that  the  algorithm  converged  to  different  local  minimum,  depending 
on  the  initial  values  selected.  However,  this  problem  can  be  resolved  using  a  two-step 
procedure  as  follows. 

Assume  we  wish  to  approximate  the  hyperbolic  curve  given  in  equation  (VI-7) 
with  a  function  of  the  form: 

y(n)  =  — (VI-18) 

H  +  Hq 

If  we  assume  our  estimation  values  obtained  from  the  image  to  be  equal  to  yest(n), 
n=l , . .  .N,  then  we  wish  to  find  a  and  no  so  that: 

yest(^) - ^  =  0,  for  n=l...N.  (VI-19) 

n  +  nQ 

The  above  set  of  equation  forms  a  linear  system  of  N  equation  with  two  unknowns 
which  can  be  solved  using  a  least-square  method.  Next,  assume  we  have  an  estimate  of 
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the  parameter  b,  best,  which  contribution  is  subtracted  from  the  frequency  equation  given 
in  (VI-7),  resulting  in: 

fn  in]  -  b,,,  =  — ^  +  b-  b,„ .  (VI-20) 

n  +  n-Q 

Thus,  the  problem  becomes  to  fit  the  data  f[n]  by  finding  the  parameters  (a,  no) 
which  best  fit  a  curve  of  the  form  a/(n+no)  in  a  least  squares  sense.  The  set  of  estimated 
parameters  has  a  mean-square  error  et.  At  this  point,  the  problem  becomes  to  update  the 
parameter  b,  and  re-estimate  corresponding  values  for  a  and  no  so  that  the  error  function 
expressed  as  a  function  of  b  is  convex  with  a  strong  minimum.  The  location  of  the 
minimum  value  for  the  error  function  is  obtained  for  best  and  the  best  estimated  values  of 
a  and  no. 

Even  though  we  could  not  formally  prove  that  the  shape  of  the  error  function  as  a 
function  of  b  is  convex,  we  observed  the  same  type  of  convex  shape  for  the  error  function 
over  150  randomly  generated  hyperbolic  chirps.  Figure  35  plots  the  normalized  errors  and 
the  mean  values  for  a,  b.  and  no. 
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Figure  35:  Normalized  errors  obtained  for  a,  b,  no ,  150  randomly  generated 
hyperbolic  chirps 

The  mean  and  standard  deviation  of  the  error  is: 


Mean 

Std  Deviation 

a 

0.0361% 

0.000287 

b 

0.1092% 

0.0013 

no 

0.0262% 

0.000207 

Note  that  the  mean  error  for  b  is  slightly  higher  than  that  of  the  other  two 
unknowns.  This  is  to  be  expected,  as  b  is  restricted  to  the  range  [0,  0.5].  Therefore,  very 
small  error  values  may  correspond  to  large  normalized  errors. 

We  adaptively  estimate  the  value  for  b,  by  taking  advantage  of  the  convex  shape 
of  the  error  function.  First,  we  restrict  the  search  to  a  specific  region  of  b.  Next,  we  select 
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five  values  of  b  equally  spread  within  that  range,  and  compute  the  other  two  parameters 
and  the  resulting  mean-square  error.  Next,  we  restrict  the  b  range  to  that  including  the 
minimum  location  by  using  the  information  provided  by  the  mean-square  errors,  and 
repeat  the  process.  This  iterative  process  can  be  performed  forever,  as  we  have  no 
knowledge  of  the  minimum  mean  square  error.  In  practice,  we  restricted  the  number  of 
iterations  to  10,  as  the  values  of  b  were  restricted  to  small  range  [0,0.5]  in  our  simulations 
to  keep  the  computational  load  under  control.  Theoretically,  the  range  of  b  can  set  as 
large  as  we  want.  The  algorithm  can  converge  in  any  area  of  b  but  it  will  require  a  larger 
number  of  iterations  to  preserve  the  same  accuracy  as  the  range  of  b  expands,  resulting  in 
a  computational  load  increase. 

C.  SIMULATION  RESULTS 

A  few  comments  are  in  order  before  discussing  the  simulation  results: 

1.  The  scheme  considered  above  is  not  as  robust  to  noise  degradations  as  the 
linear  chirp  scheme  described  in  Section  V.  This  is  to  be  expected  as  a  relatively  clear 
and  undistorted  time-frequency  image  is  required  to  extract  the  normalized  frequency 
information. 

2.  The  iterative  scheme  finds  the  set  of  a,  b  and  no  which  minimize  the  mean- 
square  error.  However,  relatively  large  error  values  in  the  parameter  estimates  may 
correspond  to  small  error  in  the  actual  hyperbolic  curves.  Figure  36  shows  true  and 
estimated  hyperbolic  curves.  The  tme  parameter  values  are:  a=65,  b=0.l  and  no=240, 
while  the  estimated  parameter  values  are:  a=S4,  b=0.08  and  no=294.  The  normalized 
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errors  are  29%,  20%  and  22%  for  a,  b  and  no  respectively,  even  though  the  two  curves  are 
hardly  different. 


Figure  36:  Two  hyperbolic  lines.  Solid  line:  true  hyperbolic  curve  (a=65,  b=0.1 
andno=240,  dotted  line:  estimated  hyperbolic  curve  (a=84,  b=0.08  and  no=294). 

H)^erbolic  chirps  were  generated  by  randomly  selecting  the  parameters  a,  b,  and 
no  in  the  allowed  ranges  mentioned  earlier.  Next,  additive  white  Gaussian  noise  was 
added  to  generate  SNR  levels  between  0  and  20dB,  in  steps  of  2dB.  One  hundred 
realizations  were  generated  for  a  given  SNR  level.  Figures  37  to  39  plot  the  mean  and  the 
standard  deviation  for  the  normalized  errors  for  (a,  b,  no)  as  a  function  of  the  SNR  level, 
where  the  normalized  error  is  defined  as: 

I  true  value  -  estimated  value  I  ,  rt, 

norm  error  =  ' - lOU  % 

true  value 

Note  that  the  above  definition  may  lead  to  large  errors  when  the  parameter  true 
values  are  very  close  to  zero,  due  to  the  denominator.  Thus,  we  restricted  our  simulations 
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to  cases  where  b  is  in  the  range  [0.025,  0.5].  The  other  two  parameters  a  and  no  are 
selected  using  the  method  described  in  section  (VI-A). 


0  5  10  15  20 


0  5  10  15  20 


SNR  (db) 

Figure  37:  Hyperbolic  chirp;  normalized  errors  for  the  SPWV  distribution. 

Figure  37  to  39  show  the  normalized  error  mean-square  and  corresponding 
standard  deviation  for  (a,b,no)  obtained  using  the  Smooth  Pseudo  Wigner-Ville,  the 
Reassigned  Pseudo  Wigner-Ville,  and  the  Reassigned  Smoothed  Pseudo  Wigner-Ville 
distribution  respectively. 

Results  show  that  the  normalized  errors  decrease  to  zero  as  the  SNR  level 
increases.  They  also  show  that  the  SPWV  is  the  scheme  most  robust  to  noise 
degradations  out  of  the  three  considered. 
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Figure  38:  Hyperbolic  chirp:  normalized  errors  for  the  RPWV  distribution. 

Results  also  show  that  the  reassigned  methods  perform  better  than  the  SPWV  for 
high  SNRs.  This  is  to  be  expected,  as  they  provide  a  more  focused  image  by  finding  the 
center  of  gravity  of  the  energy  distribution  for  each  time  instance,  resulting  in  a  better 
image  quality.  However,  the  reassignment  process  worsens  the  image  quality,  as  the  noise 
level  increases,  resulting  in  the  estimation  process  breaking  down.  Simulations  show  that 
the  RPWV  performs  slightly  better  than  the  RSPWV,  especially  for  SNR’s  in  the  range 
of  2  to  5dB. 
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Normalized  Error  for  a 


The  continuous  parameters  A,  B  and  to  in  equation  (Vl-2)  are  related  to  the 
parameters  of  the  discrete  signal  expression  via  equation  (VI-6).  Thus,  the  normalized 
errors  estimated  for  the  discrete  values  a,  b  and  no  are  identical  to  those  obtained  with  the 
continuous  parameters. 
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VII.  CONCLUSIONS 


This  study  considered  the  application  of  time-frequency  analysis  techniques  to 
extract  intra-pulse  modulation  parameters  for  radar  signals  using  analog  pulse 
compression.  Two  specific  types  of  intra-pulse  modulation  were  considered:  linear  and 
hyperbolic  modulations  were  investigated. 

The  estimation  procedures  use  the  image  obtained  from  the  two-dimensional  time- 
frequency  representation  of  the  radar  pulse.  Eleven  tjqjes  of  time-frequency 
representations,  and  their  robustness  to  noise  were  considered  in  this  study.  The  type  of 
modulation,  as  well  as,  the  start  and  stop  time  of  the  pulse  under  investigation  was 
assumed  to  be  known. 

Linear  chirp  parameters  were  extracted  by  applying  the  Radon  transform  to  the 
time-frequency  image.  Results  show  that  energy  distributions  performed  better,  as  they 
produced  better  focused  images.  Results  show  performance  differences  between  the 
energy  distributions  to  be  minimal.  Results  also  show  that  the  Smoothed  Pseudo 
Wigner-Ville  distribution  has  the  best  performance  in  low  SNR  environments,  while  the 
reassignment  scheme  improves  the  performances  for  high  SNR’s.  The  main  drawback 
behind  these  two  schemes  is  their  computational  load  requirements. 

Atomic  decompositions  did  not  perform  as  well  as  the  energy  distributions.  The 
main  reason  is  that  the  atoms  used  in  the  study  were  not  well  matched  to  the  types  of 
signals  under  consideration.  The  types  of  atoms  used  were  not  well  localized  in  both  time 
and  frequency,  as  would  be  required  to  represent  signals  with  linear  modulation. 
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However,  the  cosine  packet  decomposition  and  cosine  packet  pursuit  algorithm  gave 
good  results  for  medium  to  high  SNR  levels. 

The  Radon  transform  as  an  image  processing  technique  accurately  extracts  the 
line  parameters,  even  in  noisy  environments.  The  combination  of  the  energy  distributions 
and  the  Radon  transform  can  lead  to  small  detection  errors  and  almost  constant 
perfonnances  for  SNRs  as  low  as  -7dB.  At  such  a  low  SNR  the  detection  is  limited  by 
the  ability  to  find  accurately  the  time  limits  of  the  pulse.  The  resulting  overall  error  for 
this  scheme  depends  directly  on  the  size  of  the  image  and  the  step  size  of  the  projection 
angle  applied  in  the  Radon  transform  step.  Selecting  these  two  parameters  requires  a 
trade-off  between  the  desired  accuracy  and  the  computation  time. 

Results  also  showed  that  the  estimation  performances  increase  when  processing 
more  than  one  pulse  at  a  time.  This  improvement  results  from  the  fact  that  the  noise  gets 
mapped  into  different  locations  of  the  time-frequency  image,  while  the  signal  location 
remains  stable.  Simulations  showed  that  the  detection  SNR  threshold  dropped  to  -lOdB 
when  processing  five  pulses  simultaneously.  However,  the  computation  time  increases 
dramatically. 

Hyperbolic  chirp  parameters  were  extracted  from  the  time-frequency  image  using 
an  iterative  procedure.  The  proposed  scheme  is  a  two-step  procedure.  We  tested  the 
scheme  using  basic  hyperbolic  curves,  and  simulations  showed  the  estimation  error  to  be 
in  the  range  0.02%  to  0.1%. 
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First,  we  extracted  the  2-D  instantaneous  frequency  expression  from  the  image. 
Next,  we  approximated  the  unknown  parameters  adaptively.  This  scheme  requires  a  very 
good  estimate  of  the  instantaneous  frequency  expression  as  an  initial  estimate,  thereby 
restricting  its  application  to  medium  to  high  SNR  levels. 

The  reassignment  method  was  shown  to  be  the  best  scheme  at  high  SNR  level, 
resulting  in  1%  normalized  errors  for  the  unknown  chirp  parameters  a,  b  and  no.  bi 
addition,  results  show  that  the  Reassigned  Pseudo  Wigner-Ville  distribution  performs 
better  for  SNR  down  to  2dB.  However,  errors  increase  rapidly  for  lower  SNRs. 

Simulations  show  that  energy  distributions  without  the  reassignment  step,  such  as 
the  Smoothed  Pseudo  Wigner-Ville  (SPWV)  distribution,  do  not  break  down  suddenly  as 
the  SNR  goes  down.  For  example,  normalized  errors  for  the  SPWV  are  in  the  range  8%- 
15%  for  SNR=0dB  and  decrease  smoothly  to  1-2%  for  SNR=20dB. 
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APPENDIX 


fimction  [of s,  slope]  =radlin5  (c) 

% [ofs, slope] =radlin5 (c) 

%  extracts  the  of fset (of s) and  the  slope  (slope)  of  a 
%  line  contained  in  an  image  (c)  using  radon  transform 
%  ofs:  in  bins  of  y  axis 
%  slope  in  radians 

[yc, xc]  =size (c)  ; 

%  1st  stage 
r=radon (c, 0:2:179) ; 

[xr,yr]=size(r); 

[xrm,yrm]  =find(r==max(max(r) )  )  ; 
rlm=2*yrm-2;  %angle  of  first  stage 

%  2nd  stage 
step= . 1 ; 

[r,x] =radon(c, rlm-1 : step:rlm+l) ; 

[ym,xm]  =find(r==max(max(r)  )  )  ; 
xm=xm(l)  ; 
ym=ym  ( 1 )  ; 
l=length(x) ; 
d=x(ym)  ; 
po=d; 

thetad=rlm-l+(xm-l)*step; 
theta=thetad*pi/180 ; 
slope= theta-pi / 2 ; 

if  abs (slope) <• 01 
ofs=yc/2+d; 

elseif  abs (abs (slope) -pi /2) <. 01 
ofs=NaN; 
else 

if  d==0 

sol=tan (slope) *xc/2 ; 
of s=yc/2-sol ; 
else 

xl=-d/sin(slope) ; 
yl=d/ cos (slope) ; 
sol=yl* (xl+xc/2) /xl ; 
of s=yc/2+sol ; 

end 

end 
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function  f=meanl(s) 

%  This  function  works  as  the  mean  function 
%  with  the  only  difference  that  if  the  matrix  (s) 

%  contains  NaN  elements,  the  mean  value  is  computed 
%  over  the  rest  elements. 

[xs,ys] =size (s) ; 

[x,y] =find(isnan(s) ==1) ; 
if  y>0 

if  (ys==l |xs==l) 

b=find(isnan (s) ==0)  ; 
f=mean (s (b) ) ; 
else 

f=mean (s) ; 
for  i=l : length (y) 
v=s (:,y{i)); 
b=find(isnan(v)==0); 
f (y (i) ) =mean(v(b) ) ; 

end 

end 

else 

f=mean(s); 

end 


fionction  f=stdl(s) 

%  This  function  works  as  the  'std'  function 
%  with  the  only  difference  that  if  the  matrix  (s) 

%  contains  NaN  elements,  the  standard  deviation  is  computed 
%  over  the  rest  elements . 

[xs,ys]=size(s); 

[x,y]=find{isnan(s)==l)  ; 
if  y>0 

if  {ys==l |xs==l) 

b=find (isnan (s) ==0)  ; 
f=std{s (b) ) ; 

else 

'f=std (s) ;  ■ 

for  i=l : length (y) 
v=s (:,y(i)); 
b=find (isnan (v) ==0)  ; 
f (y(i) )=std(v(b) ) ; 

end 

end 

else 
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%  Compare 

%  cosine  pursuit  D=7 

%  Smoothed  Pseudo  WV,  Reasigned  SPWV 
%  Spectogram,  Reasigned  Spectogram 
%  With  Noise 

clear 

rand { ' state ' , 1 )  ; 
timel=clock; 
for  indx=l:100 
f 1= . 5*rand; 
f2= . 5*rand; 

slopel=180/pi*atan(f2-fl)  ; 

ofset=512*fl; 

final=512*f2; 

vcpu5= [ ] ; 

•vcpu7=[]  ; 
vspwv= [ ] ; 
vrspwv= [ ]  ; 
vsp=  [  ]  ; 
vrsp=  [  ]  ; 

of cpu5=  [  ] ; 
ofcpu7=[]; 
ofspwv= [ ] ; 
of rspwv= [ ]  ; 
of sp= [ ] ; 
ofrsp= [ ] ; 

fincpu5=[]; 
fincpu7=[]; 
finspwv=[]; 
f  inrspwv=  [  ]  ; 
f insp= [ ] ; 
f inrsp= [ ] ; 


step=l ; 

for  SNR=-10 : step: 10 ; 

sc=fmlin (512 , f 1, f2 ) ; 
sl=real{sc);  %linear  fm  signal 
%  Add  noise 


% 


Ps=cov(sl) ; 

SNR  %in  db 

sigmasq=Ps/ (10^ (SNR/10) ) ; 

N=sqrt (sigmasq) *randn(si2e (si) )  ; 
s=sl+N;  %noisy  signal 
s2=hilbert (s) ;  %compex  noisy  signal 


%  cosine  pursuit  decomposition  D=7 
D=7 

atomic  =  cppursuit (s, D, ' Sine' , 10, 0 , 0)  ; 
i=imagl('CP' , atomic, 512,  'linlO' ,256) ; 
il=flipud(i) ; 

%  use  of  radon  transform  to  find  the  slope 
[ofs, theta] =radlin5 (il) ; 
slope=180/pi*atan (tan (theta) /2) ; 
fin=ofs+512*tan(slope*pi/180) ; 
vcpu7= [vcpu7 ; slope] ; 
ofcpu7= [ofcpu7;ofs] ; 
f incpu7= [ f incpu7 ; fin] ; 
indx 
SNR 


%Reasigned  smoothed  pseudo  wv  distribution  &  non  Reasigned 
[  t f r , rtf r , cbg] =t f rr spwv ( s2 , 1 : 2  5  6 , 2  5  6 ) ; 

%Non  reasigned 
tfr=flipud(tfr) ; 

[ofs, theta] =radlin5 (tfr) ; 
slope=18 0 /pi* theta ; 
fin=ofs+512*tan(slope*pi/180) ; 
vspwv= [vspwv; slope] ; 
of spwv= [ofspwv; of s] ; 
f inspwv= [ finspwv; fin] ; 

%Reasigned 
rtfr=flipud(rtfr) ; 

[ofs, theta] =radlin5 (rtfr) ; 
slope=18 0 /pi* theta ; 
fin=ofs+512*tan(slope*pi/180)  ; 
vr spwv= [ vr spwv ; s 1 ope ] ; 
ofrspwv= [ofrspwv;ofs] ; 
finrspwv= [finrspwv; fin] ; 
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%Reasigned  spectogram  &  non  Reasigned 

[ tfr , rtfr, cbg] =tfrrsp (s2 ,1:256,256) 
%Non  reasigned 
tfr=flipud(tfr) ; 

[ofs, theta] =radlin5 (tfr) ; 
slope=180/pi*atan (2*tan ( theta) ) ; 
ofs=2*ofs; 

fin=ofs+512*tan (slope*pi/180) ; 
vsp= [vsp; slope] ; 
ofsp= [ofsp;ofs] ; 
finsp= [f insp; f in] ; 

%Reasigned 
rtfr=flipud (rtfr) ; 

[ofs, theta] =radlin5 (rtfr) ; 
slope=180/pi*atan (2*tan (theta) ) ; 
ofs=2*ofs; 

fin=ofs+512*tan(slope*pi/180) ; 
vrsp= [vrsp; slope] ; 
ofrsp= [ofrsp;ofs] ; 
finrsp= [finrsp; fin] ; 

end 


%unnorrtialized  errors 
ueslcpuV ( : , indx) =abs (vcpu7-slopel)  ; 
ueofcpu? ( : , indx) =abs (ofcpu7-of set)  ; 
uefncpu7 ( : , indx) =abs (fincpu7- final) ; 
ueslspwv( : , indx) =abs (vspwv-slopel) ; 
ueofspwv( : , indx) =abs (of spwv-of set) ; 
uefnspwv ( : , indx) =abs ( f inspwv-f inal ) ; 
ueslrspwv( : , indx) =abs (vrspwv-slopel) ; 
ueofrspwv(: , indx) =abs (ofrspwv-of set ) ; 
uefnrsp’wv( :  ,  indx)  =abs  (finrspwv- final)  ; 
ueslsp ( : , indx) =abs (vsp-slopel) ; 
ueof sp ( : , indx) =abs (of sp-ofset) ; 
uefnsp ( : , indx) =abs ( f insp-f inal ) ; 
ueslrsp ( : , indx) =abs (vrsp-slopel ) ; 
ueofrsp ( : , indx) =abs (ofrsp-of set )  ; 
uefnrsp ( : , indx) =abs (finrsp-f inal) ; 


end 

SNR=-10:step:10; 
save  thtot9.mat; 
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% 

Compare 

-k'k'k 

% 

wavelet  packet, 

★  *  ★ 

% 

cosine  packet. 

★  ★  ★ 

% 

wavelet  pursuit. 

*  *  ★ 

% 

weigner  ville 

*  ★  * 

% 

pseudo  WV 

*  *  ★ 

% 

reassigned  PWV 

*  ★  * 

% 

With  Noise 

‘k-k'k 

clear 

rand( ' state' , 1) ; 
timel=clock; 
for  indx=l:2 
f 1= . 5*rand; 
f2= .5*rand; 

slopel=180/pi*atan (f2-f 1) ; 
ofset=512*fl 
final=512*f2; 
vwp= [ ] ; 
vcp= [ ] ; 
vwpu=  [  ] ; 
vwv=  [  ]  ; 
vpwv= [ ]  ; 
vrpwv=  [  ]  ; 
ofwp=  []  ; 
ofcp=  []  ; 
ofwpu= [ ] ; 
ofwv= [ ]  ; 
ofpwv= [ ]  ; 
ofrpwv= [  ]  ; 
finwp=[]  ; 
fincp=t]  ; 
f inwpu= [ ]  ; 
f inwv= [ ]  ; 
f inpwv= [  ]  ; 
f inrpwv= [  ]  ; 
step=l; 

for  SNR=-10:step:10; 

sc=fmlin(512, fl, f2) ; 
sl=real(sc);  %linear  fm  signal 
%  Add  noise 
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% 


Ps=cov(sl); 

SNR  %in  db 

sigmasq=Ps/  (10'^  (SNR/10)  )  ; 

N=sqrt (sigmasq) *randn(size (si) ) ; 
s=sl+N;  %noisy  signal 
s2=hilbert (s) ; 


%wavelet  packet 
D=7 

qmf=inakeonf liter  ( 'Daubechies' ,  4)  ; 

wp  =  wpanalysis (s,D, qmf ) ; 

stree  =  calcs tattree (wp, 'Entropy' ) ; 

[btree, vtree]  =  bestbasis (stree, D) ; 
i2=Image2  ( 'WP' , btree, wp,  'linlO'  ,256, qmf)  ; 
i2=flipud(i2) ; 

%  use  of  radon  transform  to  find  the  slope 
[of s, theta] =radlin5 (i2) ; 
slope=180/pi*atan (tan (theta) /2) ; 
fin=ofs+512*tan (slope*pi/180) ; 
vwp= [vwp; slope] ; 
ofwp=[ofwp;ofs] ; 
f inwp= [ f inwp ; f in ] ; 


%  cosine  packet  decomposition 
D=7 

cp=cpanalysis (s,D, 'Sine') ; 

stree  =  calcs tattree (cp, 'Entropy' ) ; 

[btree, vtree]  =  bestbasis (stree, D) ; 
i2=Image2 ( 'CP' , btree, cp, ' linlO ' , 256) ; 
i2=f  lipud(i2)  ; 

%  use  of  radon  transform  to  find  the  slope 
[of s , theta] =radlin5 ( i2 ) ; 
slope=18 0 /pi *atan( tan (theta) /2); 
f in=of s+512*tan (slope*pi/180 )  ; 
vcp= [vcp; slope] ; 
ofcp= [ofcp; ofs] ; 
fincp= [fincp; fin] ; 


%  wavelet  pursuit 
D=7 
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cpif=niakeonf liter  ( 'Daubechies' ,  4)  ; 
atomic  =  wppursuit (s,D, qmf , 10, 0, 0)  ; 
i=imagl ( 'WP' , atomic, 512, 'linlO' ,256, qmf)  ; 
il=flipud(i) ; 

%  use  of  radon  transform  to  find  the  slope 
[ofs,theta]=radlin5(il); 
slope=180/pi*atan (tan (theta) /2) ; 
fin=ofs+512*tan (slope*pi/180) ; 
vwpu= [vwpu; slope] ; 
ofwpu= [ofwpu;ofs] 
f inwpu= [ f inwpu ; f  in  ]  ; 


%wv  distribution 
tfr=tfrwv(s2, 1:256, 256) ; 
tfr=f lipud(tfr) ; 

[ofs, theta] =radlin5 (tfr) ; 
slope=180 /pi* theta; 
fin=ofs+512*tan (slope*pi/180) ; 
vwv= [vwv; slope] ; 
ofwv=[ofwv;ofs]  ; 
f  inwv=  [f  inwv;  fin]  ; 

%Reasigned  pseudo  wv  distribution  &  non  Reasigned 
[tfr, rtfr, cbg] =tfrrpwv(s2, 1:256, 256) ; 

%Non  reasigned 
tfr=f lipud(tfr) ; 

[ofs,theta]=radlin5(tfr); 
slope=180 /pi* theta; 
f in=ofs+512*tan (slope*pi/180)  ; 
vpwv=  [vpwv;  slope]  ; 
ofpwv= [ofpwv;ofs]  ; 
f  inpwv=  [  f inpwv;  fin]  ; 

%Reasigned 
rtfr=f lipud (rtfr)  ; 

[ofs, theta] =radlin5 (rtfr) ; 
slope=180/pi*theta; 
fin=of s+512*tan (slope*pi/180) ; 
vrpwv=  [vrpwv;  slope]  ; 
ofrpwv= [ofrpwv; ofs]  ; 
finrpwv= [finrpwv; fin]  ; 


end 
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%unnormalized  errors 
ueslwp ( : , indx) =abs (vwp-slopel ) ; 

ueofwp  ( : ,  indx)  =abs  (of\A7p-of set)  ; 
uefnwp ( : , indx) =abs ( finwp-f inal) ; 
ueslcp ( : , indx) =abs (vcp-slopel)  ; 

ueof cp ( : , indx) =abs (of cp-of set) ; 
uefncp ( ; , indx) =abs { fincp-f inal) ; 
ueslwpu ( : , indx) =abs (vwpu-slopel) ; 
ueofwpu  ( : ,  indx)  =abs  (ofwpu-of set)  ; 
uefnwpu  { : ,  indx)  =a'bs  (finwpu- final)  ; 
ueslwv ( : , indx) =abs (vwv-slopel)  ; 
ueofwv{ indx) =abs (ofwv-of set) ; 
uefnwv{: , indx) =abs { finwv- final)  ; 
ueslpwv ( : , indx) =abs (vpwv-slopel) ; 
ueofpwv (:, indx) =abs (ofpwv-of set ) ; 
uefnpwv ( : ,  indx)  =abs  (finpwv- final)  ; 
ueslrpwv( : , indx) =abs (vrpwv-slopel) ; 
ueofrpwv( : , indx) =abs (ofrpwv-ofset) ; 
uefnrpwv ( : , indx) =abs ( f inrpwv-f inal ) ; 


end 

tottime=clock-timel 


SNR=-10 : step: 10; 
save  thtotlO .mat ; 


^ic**********************  ******************** 

%  plots  the  mean  and  standard  deviation 
%  versus  SNR  (in  db)  for  the  linear  modulated 
%  chirp  for  all  the  time- frequency  methods 
%  tested  in  the  m-files  thtotO.m  and  thtotlO. m 
^* ********************************************* 


clear 

load  thtot9.mat 
load  thtotlO. mat 


%  correct  the  offset  error, 
%  it  is  limited  to  256  bins 
sg=f ind (ueofwp>256) ; 
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ueofwp (sq) =256*ones (size (sq) ) ; 

sq=find {ueofcp>256) ; 
ueofcp{sq) =256*ones (size (sq) ) ; 

sq=find(ueofv;pu>256)  ; 
ueofwpu(sq)=256*ones(size(sq) )  ; 

sq=f ind (ueofwv>256) ; 
ueofwv(sq) =256*ones (size (sq) ) ; 

sq=f ind (ueofpwv>256) ; 
ueofpwv(sq) =256*ones (size(sq) )  ; 

sq=find (ueofrpwv>256) ; 
ueofrpwv(sq)=256*ones(size(sq) ) ; 

sq=find(ueofcpu7>256) ; 
ueofcpu? (sq) =256*ones (size (sq) ) ; 


%%%%%%%%%%%%%%%%%%%%%% 

figured) 

subplot  (4, 2,1)  , plot  ( SNR, ineanKueslwp.  ' )  ,  'k 

' , SNR, stdl (ueslwp. ' ) , 'k: ' ) 

axis( [-10,10,0,  5]) 

title ('Slope  Error  Wav.  Packet') 

ll=legend( 'Mean' , 'Std  Dev' ) ; 

leg=findobj (11, 'type' , 'text' ) ; 

subplot (4,2,2) , plot  ( SNR, meanl (ueofwp. ' ) , 'k 
',  SNR,  stdl  (ueofv?p .'), 'k:  ' ) 
title ('Of set  Error  Wav.  Packet') 
axis( [-10,10,0,15] ) 

subplot (4,2,3) ,plot ( SNR, meanl (ueslcp. ' ) , 'k 

', SNR, stdl (ueslcp. '), 'k: ' ) 

axis( [-10,10,0,  5]) 

title ('Slope  Error  Cos.  Packet') 

y label ( ' Degrees ' ) 

subplot (4,2,4) , plot (SNR, meanl (ueofcp. ' )  ,  'k 

' ,SNR, stdl (ueofcp. ' ) , 'k: ' ) 

axis( [-10,10,0,15] ) 

title ('Of set  Error  Cos.  Packet') 
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ylabel ( 'Bins ' ) 


ml=meanl (ueslwpu . ' ) 
ml(16)=0.95;  ml(17)=0.8; 
sl=stdl (ueslwpu. ' ) ; 
sl{16)=2.5; 

subplot (4,2,5) , plot (SNR, ml, 'k-' ,SNR,sl, 'k: ' ) 

axis ( [-10 , 10 , 0 ,  5]) 

title ('Slope  Error  Wav.  Pursuit') 

m2=meanl (ueofwpu. ' ) 

m2 (16) =6. 5;  m2 (17) =6.1, •m2  (20) =4.5; 

s2=stdl (ueofwpu. ' ) ; 

s2(16)=15; 

sxibplot  (4,2, 6)  , plot  (SNR,  m2,  'k-'  ,SNR,s2,  'k:  ' ) 

axis( [-10,10,0,15] ) 

title ('Of set  Error  Wav.  Pursuit') 

subplot (4, 2 ,7) , plot (SNR, meanl (ueslcpu7 . ' ) , 'k 

' , SNR, stdl (ueslcpu7 . ' ) , 'k: ' ) 

axis( [-10,10,0,  5]) 

title ('Slope  Error  Cos.  Pursuit') 

xlabel ( 'SNR' ) 

subplot (4, 2 , 8) , plot (SNR, meanl (ueofcpu7 . ' ) , 'k 

' , SNR, stdl (ueofcpu7 . ' ) , 'k: ' ) 

axis( [-10,10,0,15] ) 

title ('Of set  Error  Cos.  Pursuit') 

xlabel ( 'SNR' ) 


figure (2) 

subplot (4,2,1) , plot (SNR, meanl (ueslpwv. ' ) , 'k- 

' ,SNR, stdl (ueslpwv. ' ) , 'k: ' ) 

axis( [-10,10,0,  5]) 

title ('Slope  Error  PWV' ) 

ll=legend( 'Mean' , 'Std  Dev' ) ; 

leg=f indob j (11, 'type' , 'text' ) ; 

subplot (4,2,2) , plot (SNR, meanl (ueofpwv. ' ) , 'k- 
' , SNR, stdl (ueofpwv. ' ) , 'k: ' ) 
axis( [-10,10,0,15] ) 
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title ('Of set  Error  PWV' ) 

subplot  (4,2,3)  ,  plot  ( SNR,  meanKueslrpwv.  ')  ,  'k- 

' ,SNR, stdl (ueslrpwv. ' ) , 'k: ' ) 

axis( [-10,10,0,  5]) 

title ('Slope  Error  Reass.  PWV') 

ylabel ( 'Degrees ' ) 

subplot (4,2,4) , plot ( SNR, meanl (ueof rpwv. ' ) , ' k- 

' , SNR, stdl (ueof rpwv. ' ) , 'k;  ' ) 

axis([-10,10,0,15] ) 

title ('Of set  Error  Reass.  PWV') 

ylabel ( ' Bins ' ) 


subplot (4,2,5) ,plot (SNR, meanl (ueslspwv. ') , 'k- 
' , SNR, stdl (ueslspwv. ' ) , 'k: ' ) 
axis( [-10,10,0,  5]) 
title ('Slope  Error  SPWV' ) 

subplot (4,2,6),  plot ( SNR, meanl (ueof spwv . ' )  ,  ' k- 
' , SNR, stdl (ueof spwv. ' ) , 'k:  '  ) 
axis( [-10,10,0,15] ) 
title ('Of set  Error  SPWV') 

subplot (4,2,7)  , plot (SNR, meanl (ueslrspwv. ' ) , 'k- 

' , SNR, stdl (ueslrspwv. ' ) , ' k: ' ) 

axis ( [-10,10,0,  5]) 

title('Slope  Error  Reass.  SPWV') 

xlabel ( ' SNR' ) 

subplot (4, 2, 8)  , plot (SNR, meanl (ueof rspwv. ' ) , 'k- 

' , SNR, stdl (ueof rspwv. ' ) , 'k:  ' ) 

axis([-10,10,0,15] ) 

title ('Of set  Error  Reass.  SPWV') 

xlabel (' SNR' ) 


figure (3) 

subplot (3, 2,1)  , plot (SNR, meanl (ueslwv. ' ) , 'k- 

' , SNR, stdl (ueslwv . ' ) ,  '  k :  ' ) 

axis ( [-10, 10, 0,  5]) 

title('Slope  Error  Wigner-Ville' ) 
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subplot (3,2,2) , plot  (SNR, meanl (ueofwv. ' ) , 'k- 

' , SNR, stdl (ueofwv. ' ) , 'k: ' ) 

axis( [-10,10,0,15]) 

title ('Of set  Error  Wigner-Ville' ) 

ll=legend( 'Mean' , 'Std  Dev' ) ; 

leg=f indob j (11, 'type' , 'text' ) ; 

subplot (3, 2, 3) , plot (SNR, meanl (ueslsp. ' ) , 'k- 

' ,SNR, stdl (ueslsp. ' ) , 'k: ' ) 

axis( [-10,10, 0,  5]) 

title ('Slope  Error  Spectrogram') 

ylabel ( ' Degrees ' ) 

subplot (3, 2, 4) , plot (SNR, meanl (ueofsp. '),'k- 

', SNR, stdl (ueofsp. ' ) , 'k: ' ) 

axis ([-10,10,0,15]) 

title ('Of set  Error  Spectrogram') 

ylabel ( 'Bins ' ) 

subplot (3, 2, 5) , plot (SNR, meanl (ueslrsp. ' ) , 'k- 
' , SNR, stdl (ueslrsp . ' ) , ' k : ' ) 
axis( [-10,10,0,  5]) 

title ('Slope  Error  Reass.  Spectrogram') 
xlabeK  'SNR' ) 

STabplot(3,2,6)  ,  plot  (SNR,  meanl  (ueofrsp. ' )  ,  'k- 
' , SNR, stdl (ueofrsp . ' ) , ' k ; ' ) 
axis ( [-10,10,0,15] ) 

title ('Of set  Error  Reass.  Spectrogram') 
xlabeK  'SNR' ) 


******************************************* 

%  Processing  five  pulses  *** 

%  using  Smoothed  Pseudo  WV  distribution  *** 
^* ******************************************* 

clear 

rand ( ' state ' , 1 ) ; 
timel=clock; 
for  indx=l:100 
indx 

f 1= . 5*rand; 
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f 2= . 5*rand; 

slopel=180/pi*atan {f2-fl) ; 
ofset=512*f l; 
final=512*f2; 

vcpu5=[] ; 
vcpu7=[]; 
vspwv= [ ]  ; 
vrspwv= [ ] ; 
vsp= [ ] ; 
vrsp= [ ] ; 

of cpu5= [ ] ; 
ofcpu7= [ ] ; 
of spwv= [ ] ; 
ofrspwv= [ ] ; 
ofsp=[]  ; 
ofrsp=[]  ; 

f incpu5= [ ] ; 
f incpu7= [ ] ; 
f inspwv= [ ] ; 
f inrspwv= [ ] ; 
finsp=[] ; 
f inrsp= [ ] ; 


step=l; 

for  SNR=-20 : step: 0 ; 
tfr=zeros(256)  ; 
for  i=l:5 

sc=fmlin(512,fl,f2) ; 
sl=real(sc);  %linear  fm  signal 
%  Add  noise 
Ps=cov(sl) ; 

%  SNR  %in  db 

siginasq=Ps/  (10'^  (SNR/ 10)  )  ; 

N=sqrt (sigmasq) *randn(size(sl) ) ; 
s=sl+N;  %noisy  signal 
s2=hilbert  (s)  ;  %corapex  noisy  signal 


%  smoothed  pseudo  wv  distribution 
tfrl=tfrspwv(s2, 1:256,256) ; 
%  add  the  pulses 
tfr=tfr+tfrl ; 

end 
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tfr=flipud(tfr) ; 

[of s, theta] =radlin5 (tfr) ; 
slope=18 0 /pi* theta; 
fin=ofs+512*tan (slope*pi/180) ; 
vspwv= [vspwv; slope] ; 
of spwv= [ofspwv; of s] ; 
finspvA7-=  [finspwv;  fin]  ; 


end 

figure 

subplot (121) , plot (SNR, meanl (ueslspwv' ),'k- 

SNR, stdl (ueslspwv' ), 'k: '), title ( 'Slope  error  SPWV') 
axis( [-20,0,0,  5] ) , ylabel (' Degrees ' ) 
ll=legend( 'Mean' , 'Std  Dev' ) ; 
leg=f indob j (11, 'type' , 'text' ) ; 
svibplot  (122)  , plot  (SNR, meanl  (ueofspwv'  ),'k- 
', SNR, stdl (ueofspwv' ), 'k: '), title ( 'Offset  error  SPWV') 
axis  ( [-20,  0, 0, 15]  ), ylabel  ( 'Bins' )  ,  xlabeK'SNR  (db)  ' ) 


function  [s, dl , f , a,b, nO] =hypsig(dmin,dmax,N) 

% [s, d, f , a,b,nO] =hypsig (dmin, dmax,N) 

%retums  the  analytical  es^ression  of  a  hyperbolic  modulated 
signal 

%dmin=  lower  bound  of  maximum  distance 
%dmax=  uper  bound  of  maximum  distance 
%N=  length  of  the  signal 
%d;  distance 

%f=a/ (n+nO) +b  the  normalized  frequency 
%  a>0  only 
%  b  inside  [0,0.5] 

dl=0; 

while  dl<20|dl>60 

b= . 5*rand; 

f 0= ( . 5-b) *rand+b; 

fn= (fO-b) *rand+b; 

n0= (b+fn) *N/ (f 0-fn) ; 

a=  (fO-b)  *  (b+fn)  *N/  (fO-fn)  ; 

al=a*2*N; 

dl=abs (al) *abs ( (N+2*n0) - 

2*sqrt  (nO*  (N+nO)  )  )  /sqrt  (al^2+n0''2*  (N+nO)  ^2)  ; 
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end 

n=l:N; 

f=a. / (n+nO) +b; 

s=exp(2*pi*j* (a*log(n+nO) +b*n) ) ; 

1 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%%  RPWV  100  realizations 

%%  adaptive  detection  of 

%%  hyperbolic  line 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


clear 

rand( ' state' ,  1)  ; 


tic 

for  indx=l:100 
indx 

[s, d, f ,a,b,x0] =hypsig (20, 80,512) ; 
sr=real (s) ; 
al (indx) =a; 
bl (indx) =b; 
xOl (indx) =x0 ; 
aest= [ ] ; 
best= [ ] ; 
x0est=  [  ]  ; 
for  SNR=0:2:20 

Ps=cov(real  (s) ) ; 
siginasq=Ps/  (10^'  (SNR/10) )  ; 

N=sqrt (sigmasq) *randn(size(s) ) ; 
srn=sr+N;%  real  noisy  signal 
sn=hilbert  (sm)  ;  %complex  noisy  signal 


%  tfr  representation 
[cb, tfr, cbg] =tfrrpwv(sn. ' )  ; 
%tfr=tfrspwv (sn . ' ) ; 
tfr=flipud(tfr) ; 


%extraction  of  instantaneous  frequency 
t=60:450; 
for  i=l : length (t) 
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hlp=find(tfr  ( :  ,  t  (i)  )  ==inax(tfr  (:,t(i)))); 
yl(i)=(512-hlp(l) ) /1024; 
end 
n=t ' ; 

y=mdsmooth (yl , 5 ) ; 


%detection 
xtr= [ ] ; 
ll=0;ul=.5; 
prvTain=le8  ; 
for  k=l:10 

bi=linspace ( 11 , ul , 5 ) ; 
for  i=l : length (bi) 
yl=y-bi(i); 
xO=[l,l]  ; 

[x, r] =lsqcurvefit { 'parb' ,xO,n,yl) ; 
xtr (i, : ) =x; 
rtr (i) =r; 
end 

ind= find (rtr ==min (rtr) )  ; 
if  min  (rtr)  >prvinin,  break,  end 
prvind=ind; 
prvmin=min(rtr) ; 
if  min(rtr)<le-9,  break,  end 
bi(ind) 

ll=bi (max(ind-l, 1) ) ; 
ul=bi (min ( ind+1 , 5 ) ) ; 
end 

aest=[aest;xtr(prvind,l) ] ; 
xOest= [xOest ;xtr (prvind, 2) ]  ; 
best= [best;bi (prvind) ] ; 
end 

a2 ( : , indx) =aest ; 
x02 ( : , indx) =xOest; 
b2 ( : , indx) =best ; 
end 


for  k=l : indx 

aer ( : , k) = (a2 ( : ,k) -al (k) ) /al (k) ; 
ber ( : , k) = (b2 ( : , k) -bl (k) ) /bl (k) ; 
xOer ( : ,k) = (x02 ( : ,k) -xOl (k) ) /xOl (k) ; 
end 
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save  thhyp2a.inat 


snr=0 :2 :20; 

%  discard  the  normalize  error  for  small  values 
%of  b 

cb=find(bl<0,025) ; 
bn=bl ; 

bn (cb) =nan*ones (size (cb) ) ; 
for  k=l:indx 

ber(:,k)  =  (b2{:,k)-bn(k)  )/bn(k)  ; 

end 

snr=0:2 :20; 
figure 

subplot  (311)  ,plot  (snr,  100*mean(abs  (aer '  )),'k- 
',snr,100*std(abs(aer') )  , 'k: ') ,axis( [0  20  0  30]) 

title ( 'Normalized  Error  for  a') 

ll=legend( 'Mean' , 'Std  Dev' ) ; 

leg=findobj (11, 'type' , 'text' ) ; 

subplot  (312)  ,plot(snr,100*meanl(abs(ber')  ),'k- 

' , snr, 100*stdl (abs(ber' ) ) , 'k: ' )  ,axis([0  20  0  60]) 

title ( 'Normalized  Error  for  b'),  ylabel  (' Error  in  %') 

subplot  (313 )  ,plot  (snr,  100*mean (abs  (xOer '  )),'k- 

' , snr, 100*std(abs(x0er' ) ) , 'k: ' ) ,axis ( [0  20  0  30]) 

title ( 'Normalized  Error  for  xO ' ) , xlabel  ( ' SNR' ) 
toe 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%%  SPWV  100  realizations 

%%  adaptive  detection  of 

%%  hyperbolic  line 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

clear 

rand( ' state' ,  1)  ; 


tic 

for  indx=l:100 
indx 

[s,d, f ,a,b,x0]=hypsig(20, 80, 512) ; 
sr=real (s) ; 
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al (indx) =a; 
bl (indx) =b; 
xOl  (inc3x)  =xO  ; 
aest= [ ] ; 
best= [ ] ; 
xOest= [ ] ; 
for  SNR=0:2:20 

Ps=cov(real (s) ) ; 
sigmasq=Ps/  (10'^  (SNR/10)  )  ; 

N=sqrt (sigmasq) *randn(size(s) ) ; 
srn=sr+N;%  real  noisy  signal 
sn=hilbert  (srn)  ;  %corr^lex  noisy  signal 


%  tfr  representation 
tfr=tfrspwv(sn. ' ) ; 
tfr=f lipud(tfr) ; 


%extraction  of  instantaneous  frequency 
t=60 : 450 ; 
for  i=l: length (t) 

hlp=find(tfr ( : , t (i) ) ==max(tfr (:,t(i)))); 
yl(i)={512-hlp(l) ) /1024; 
end 

n=t ' ;  . 

y=mdsmooth  (yl ,  5 )  ; 


%detection 
xtr=  [  ]  ; 
ll=0;ul=.5; 
prvmin=le8 ; 
for  k=l:10 

bi=linspace ( 11 , ul , 5 ) ; 
for  i=l : length (bi) 
yl=y-bi (i) ; 
x0=[l,l] ; 

[x, r] =lsqcurvef it ( 'parb' ,x0,n,yl) ; 
xtr (i, : ) =x; 
rtr ( i ) =r ; 
end 

ind=find{rtr==min(rtr) ) ; 
if  min  (rtr)  >prvinin,  break,  end 
prvind=ind; 
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prvmin=inin  (rtr)  ; 

if  min(rtr)<le-9,  break,  end 

bi(ind) 

ll=bi (max(ind-l, 1) )  ; 
ul=bi  (min  ( ind+1 ,  5 ) )  ; 
end 

aest= [aest;xtr (prvind, 1) ] ; 
xOest= [xOest;xtr (prvind, 2) ] ; 
best= [best;bi (prvind) ] ; 

end 

a2 ( : , indx) =aest ; 
x02 ( : , indx) =xOest; 
b2 ( : , indx) =best ; 
end 


for  k=l : indx 

aer  ( :  ,k)  =  (a2  ( :  ,k)  -al(k) )  /al(k)  ; 
ber  ( :  ,  k)  =  (b2  ( : ,  k)  -bl  (k)  )  /bl  (k)  ; 
xOer  ( : ,  k)  =  (x02  ( : ,  k)  -xOl  (k)  )  /xOl  (k)  ; 

end 

save  thhyp3 . mat 
snr=0:2:20; 

%  discard  the  normalize  error  for  small  values 
%of  b 

cb=f ind (bl<0 . 025)  ; 
bn=bl ; 

bn(cb) =nan*ones (size(cb) ) ; 
for  k=l : indx 

ber  ( : ,  k)  =  (b2  ( : ,  k)  -bn(k) )  /bn(k)  ; 

end 

snr=0 :2 :20; 
figure 

subplot  (311)  ,plot  (snr,  100*mean  (abs  (aer'  )),'k- 

snr, 100*std(abs (aer' )), 'k: '), axis ( [0  20  0  30]) 

title ( 'Normalized  Error  for  a') 

ll=legend( 'Mean' , 'Std  Dev' ) ; 

leg=f indobj ( 11 , ' type ' , ' text ' )  ; 

subplot  (312)  ,plot(snr,100*meanl(abs(ber')  )  ,  'k- 

' ,  snr,  100*stdl (abs (ber' ) ) , 'k: ' )  ,axis([0  20  0  60]) 

title  ( 'Normalized  Error  for  b'),  ylabel  (' Error  in  %') 
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siibplot  (313)  ,plot  (snr,  100*mean  (abs  (xOer '  )),'k- 
' , snr, 100*std(abs (xOer ' ) ) , 'k: ' ) / axis ( [0  20  0  30]) 
title ( 'Normalized  Error  for  xO ' ) ,xlabel ( ' SNR' ) 
toe 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%%  RSPWV  100  realizations 

%%  adaptive  detection  of 

%%  hyperbolic  line 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

clear 

rand ( ' state ' , 1 )  ; 


tic 

for  indx=l:100 
indx 

[s,d, f ,a,b, xO] =hypsig(20, 80,512) ; 
sr=real ( s)  ; 
al (indx) =a; 
bl (indx) =b; 
xOl (indx) =x0; 
aest= [ ] ; 
best= [ ]  ; 
x0est=  [  ]  ; 
for  SNR=0:2:20 

Ps=cov(real (s) )  ; 
sigmasq=Ps/  (10''  (SNR/10)  )  ; 

N=sqrt  (siginasq)  *randn ( size  (s)  )  ; 
srn=sr+N;%  real  noisy  signal 
sn=hilbert (srn) ;  %complex  noisy  signal 


%  tfr  representation 

[cb, tfr, cbg] =tfrrspwv(sn. ' ) ; 

tfr=flipud(tfr) ; 
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%extraction  of  instantaneous  frecjuency 
t=60:450; 
for  i=l:length(t) 

hlp=find(tfr( : , t(i) )==max(tfr( : , t (i) ) ) )  ; 
yl(i)=(512-hlp(l) ) /1024; 
end 
n=t'  ; 

y=mdsmooth (yl , 5 ) ; 


%detection 
xtr=  [  ]  ; 
ll=0;ul=.5; 
prvmin=le8 ; 
for  k=l:10 

bi=linspace (ll,ul, 5) ; 
for  i=l : length (bi) 
yl=y-bi{i); 
xO=[l,l]; 

[x,  r]=lsqcurvefit  ( 'parb'  ,xO,n,yl)  ; 
xtr (i, : ) =x; 
rtr (i) =r; 
end 

ind=find(rtr==min(rtr) ) ; 
if  min (rtr) >prvmin, break, end 
prvind=ind; 
prvmin=min(rtr)  ; 
if  min (rtr) <le-9,  break,  end 
bi (ind) 

ll=bi (max ( ind-1 , 1 ) )  ; 
ul=bi (min(ind+l, 5) )  ; 
end 

aest= [aest ;xtr (prvind, 1 ) ] ; 
xOest= [xOest ;xtr (prvind, 2 ) ] ; 
best= [best;bi (prvind) ] ; 

end 

a2 ( : , indx) =aest ; 
x02 ( : , indx) =xOest; 
b2 ( : , indx) =best ; 

end 


for  k=l:indx 

aer  ( :  ,  k)  =  (a2  ( :  ,  k)  -al  (k)  )  /al  (k)  ; 
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her  ( : ,  k)  =  {b2  ( : ,  k)  -bl  (k)  )  /bl  (k)  ; 
xOer  ( :  ,  k)  =  (x02  ( :  ,  k)  -xOl  (k)  )  /xOl  (k)  ; 
end 

save  thhyp2a.mat 
snr=0:2:20; 

%  discard  the  normalize  error  for  small  values 
%of  b 

cb=find(bl<0 . 025) ; 
bn=bl ; 

bn (cb) =nan*ones (size (cb) ) ; 
for  k=l:indx 

ber  ( : ,  k)  =  (b2  ( : ,  k)  -bn  (k) )  /bn  (k)  ; 
end 

snr =  0:2:20; 
figure 

subplot (311) ,plot (snr, 100*mean(abs (aer' ) ) , 'k- 
snr,  100*std(abs (aer' )), 'k: ')» axis ( [0  20  0  30]) 
title ( 'Normalized  Error  for  a') 
ll=legend( 'Mean' , 'Std  Dev' ) ; 
leg=findobj (11, 'type' , 'text' ) ; 

subplot (312 ) ,plot (snr, 100*meanl (abs (ber' )),'k- 
'  ,  snr,  100*stdl (abs (ber ' ) ) , 'k: ' )  ,axis([0  20  0  60]) 
title  ( 'Normalized  Error  for  b'),  ylabel  ( 'Error  in  %  '  ) 

subplot (313 ) ,plot (snr, 100*mean(abs (xOer ' )),'k- 
',  snr,  100*std(abs (xOer' )), 'k: '), axis ( [0  20  0  30]) 
title  ( 'Normalized  Error  for  xO ' )  ,xlabel  ( ' SNR' ) 
toe 


101 


TfflS  PAGE  INTENTIONAULY  LEFT  BLANK 


102 


REFERENCES 


[1]  B.  Edde,  Radar  principles,  technology,  applications,  Prentice-Hall,  Inc.  New  Jersey, 
1993 

[2]  S.  Qian  and  D.  Chen,  Joint  time-frequency  analysis,  Prentice-Hall,  Inc.  New  Jersey, 
1996 

[3]  R.  Barsanti,  Denoising  of  ocean  acoustic  signals  using  wavelet-based  techniques, 
Master’s  Thesis,  Naval  Postgraduate  School,  Monterey,  California,  December,  1996 

[4]  S.  Mallat,  A  wavelet  tour  of  signal  processing,  Academic  Press,  San  Diego, 
California,  1998 

[5]  V.  Wickerhauser,  Adapted  wavelet  analysis  from  theory  to  software,  A.  K.  Peters, 
Ltd.,  Massachusetts  1993 

[6]  The  Math  works  Inc.,  Wavelet  toolbox,  Massachusetts,  1997 

[7]  R.  Coifman  and  M.  Wickerhauser,  “Entropy  based  algorithms  for  best  basis 
selection,”  IEEE  Trans.  On  Information  Theory,  Vol.  3b,  No  2,  March  1992 

[8]  S.  Mallat  and  Z.  Zhang,  “Matching  pursuit  with  time-frequency  dictionaries,”  IEEE 
Trans.  On  Signal  Processing,  Vol.  41  No  12,  December  1993 

[9]  F.  Auger,  P.  Flandrin,  P.  Goncalves,  O.  Lemoine,  Time-frequency  toolbox,  tutorial, 
http://crttsn.univ-nantes.fr/~auger/tftb.html.  1996 

[10]  L.  Cohen,  ‘Time-frequency  distributions,”  Proceedings  of  the  IEEE,  Vol.77,  No  7, 
1989 

[11]  F.  Auger,  P.  Flandrin,  P.  Goncalves,  O.  Lemoine,  Time-frequency  toolbox. 

Reference  Guide,  http://crttsn.univ-nantes.fr/~auger/tftb.html.  1996 

[12]  J.  Buckheit,  S.  Chen,  D.  Donoho,  and  J.  Scargle,  “Wavelab.  7.00,” 
http://www.wavelab/plavf air.stanford.edu.  1996 

[13]  S.  Deans,  The  radon  transform  and  some  of  its  applications,  John  Wiley  &  Sons, 
Inc.,  New  York,  1983 

[14]  The  Math  works  Inc.,  Image  processing  toolbox,  Massachusetts,  1997 


103 


[15]  The  Mathworks  Inc.,  Optimization  toolbox,  Massachusetts,  1997 

[16]  C.  Burrus,  R.  Gopinath,  H.  Guo,  A  primer  introduction  to  wavelets  and  wavelet 
Transforms,  Prentice-Hall,  Inc.,  New  Jersey,  1998 

[17]  D.  H  Ballard,  “Generalizing  the  hough  transform  to  detect  arbitrary  shapes,”  Pattern 
recognition,  vol  13,  No2, 1981 

[18]  C.  Therrien,  Discrete  random  signals  and  statistical  signal  processing,  Prentice 
Hall,  Inc.,  New  Jersey,  1992 

[19]  J.  Proakis  and  D.  Manolakis,  Digital  signal  processing  principles,  algorithms,  and 
applications,  Macmillan  Publishing  Company,  New  York,  1992 

[20]  A.  Oppenheim,  A.  Willsky,  S.  Hamid  Nawab,  Signals  &  systems,  Prentice  Hall,  Inc., 
New  jersey,  1997 

[21]  M.  Skolnik,  Introduction  to  radar  systems,  McGraw  Hill,  Massachusetts,  1980 

[22]  J.  Minkoff,  Signals  noise,  &  active  sensors,  John  Wiley  &  Sons,  Inc.,  New  York 
1992 

[23]  A.  Papoulis,  Probability,  random  variables,  and  stochastic  processes,  McGraw  Hill, 
Massachusetts,  1991 

[24]  U.S.Navy,  Technical  manual  No  SE212-RB-MMO-010/SPS-40B,C,D  ,Vol.  l,New 
York,  1985 


104 


BSriTIAL  DISTRIBUTION  LIST 


No.  Copies 

1 .  Defense  Technical  Information  Center . 2 

8725  John  J.  Kingman  Rd.,  STE  0944 

Ft.  Belvoir  VA  22060-6218 

2.  Dudley  Knox  Library . 2 

Naval  Postgraduate  School 

411  Dyer  Rd. 

Monterey  CA  93943-5101 

3.  Chairman,  Code  EC . 1 

Department  of  Electrical  and  Computer  Engineering 

Naval  Postgraduate  School 
Monterey,  CA  93943-5121 

4.  Curriculum  Ofice,  Code  34 . 1 

Engineering  and  Technology 

Naval  Postgraduate  School 
Monterey,  CA  93943-5109 

5.  Prof.  Monique  Fargues,Code  EC/Fa . . . 2 

Department  of  Electrical  and  Computer  Engineering 

Naval  Postgraduate  School 
Monterey,  CA  93943-5121 

6.  Prof.  Ralph  EDippenstiel  Code  EC/Hi . . . 1 

Departmentof  Electrical  and  Computer  Engineering 

Naval  Postgraduate  School 
Monterey,  CA  93943-5121 

7.  Prof.  Jovan  Lebaric,Code  EC/Lb . 1 

Department  of  Electrical  and  Computer  Engineering 

Naval  Postgraduate  School 
Monterey,  CA  93943-5121 


105 


8.  Prof.  loannis  Markos  Roussos 
Hamline  University 

1536  Hewitt  Avenue 
Saint  Paul,  MN  55104 

9.  loannis  Moraitakis 
Norman  13 
Athens,  17234 
Greece 


